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ZUSAMMENFASSUNG 

Obwohl Staub nur etwa ein Prozent der Masse einer protoplanetaren Scheibe aus- 
macht, ist er doch ein wichtiger Bestandteil von chemischen Modellen, Modellen 
zur Planetenentstehung oder in der Modellierung von Strahlungstransport und 
Beobachtungen. Der anfangliche Wachstumsprozess von sub-/im groBen Staub- 
partikeln bis bin zu Planetesimalen, sowie der radiale Transport von Staub in den 
Gasscheiben um junge Sterne ist das Thema dieser Arbeit. Radiale Drift, Sedi- 
mentation, turbulenter Transport sowie Koagulation, Fragmentation und Erosion 
bestimmen die zeitliche Entwicklung von zirkumstellarem Staub. 

Wir gehen dieses Problem von drei verscliiedenen Riclitungen an: analytische 
Berechnungen, numerische Simulationen und Vergleich zu Beobachtungen. Wir 
besclireiben die physikalischen und numerisclien Konzepte, mit denen unser Mod- 
ell die Entwicklung des Staubes liber Millionen von Jaliren in einer zeitabhangi- 
gen, viskosen Gassclieibe simuliert. Wir vergleichen die simulierten StaubgroBen- 
verteilungen mit unseren analytischen Vorhersagen und leiten daraus ein ein- 
faclies Rezept zur Bestimmung von stationaren GroBenverteilungen ab. Mit dem 
vorgestellten Modell ist es uns moglich zu zeigen, dass Fragmentation erklaren 
kann, warum zirkumstellare Scheiben fiir melirere Millionen Jahre reicli an Staub 
sein konnen. Des Weiteren befassen wir uns mit dem Problem, vor das uns 
Beobachtungen stellen, namlich dass Staubpartikel von der GroBe einiger Millime- 
ter in groBen Abstanden von ihrem Zentralstern nachgewiesen wurden. Unter der 
Annahme, dass radiale Drift ineffektiv ist, konnen wir einige beobachtete spektrale 
Indices und Fliisse im mm-Wellenlangenbereich reproduzieren. Lichtschwachere 
Quellen deuten darauf hin, dass das Staub-zu-Gas Verhaltnis oder die Opazitaten 
geringer sind als weithin angenommen. 



Abstract 

Dust constitutes only about one percent of the mass of circumstellar disks, yet 
it is of crucial importance for the modeling of planet formation, disk chemistry, 
radiative transfer and observations. The initial growth of dust from sub-/im sized 
grains to planetesimals and also the radial transport of dust in disks around young 
stars is the topic of this thesis. Circumstellar dust is subject to radial drift, vertical 
settling, turbulent mixing, collisional growth, fragmentation and erosion. 

We approach this subject from three directions: analytical calculations, nu- 
merical simulations, and comparison to observations. We describe the physical 
and numerical concepts that go into a model which is able to simulate the radial 
and size evolution of dust in a gas disk which is viscously evolving over several 
million years. The resulting dust size distributions are compared to our analytical 
predictions and a simple recipe for obtaining steady-state dust size distributions is 
derived. With the numerical model at hand, we show that grain fragmentation can 
explain the fact that circumstellar disks are observed to be dust-rich for several 
million years. Finally, we investigate the challenges that observations present to 
the theory of grain evolution, namely that grains of millimeter sizes are observed 
at large distances from the star. We have found that under the assumption that 
radial drift is ineffective, we can reproduce some of the observed spectral indices 
and fluxes. Fainter objects point towards a reduced dust-to-gas ratio or lower dust 
opacities. 



The light dove, cleaving in free flight the thin air, 
whose resistance it feels, might imagine, that her 
movements would be far more free and rapid in 
airless space. 

Critique of pure reason (1787), Introduction, III. 

- Immanuel Kant 



Contents 



Contents 11 

1 Introduction 15 

1.1 Historical perspective 15 

1.2 Theoretical concepts 18 

1.2.1 Structure and evolution of circumstellar disks 19 

1.2.2 The dynamics of dust 23 

1.2.3 From dust to planets 26 

1.3 Observational constraints 28 

1.3.1 Disk observations 28 

1.3.2 Meteorites 31 

1.3.3 Exoplanets 33 

1.4 Summary 34 

1.5 About this thesis 35 

2 Evolution of gas and dust in protoplanetary disks 39 

2.1 Introduction 40 

2.2 Model 41 

2.2.1 Evolution of gas surface density 42 

2.2.2 Radial evolution of dust 44 

2.2.3 Temperature and vertical structure 47 

2.2.4 Opacity 49 

2.2.5 Initial infall phase: cloud collapse 49 

2.2.6 Grain growth and fragmentation 51 

2.3 Resuhs 58 

2.3.1 Viscous evolution of the gas disk 58 

2.3.2 Fiducial model without fragmentation 59 

2.3.3 Fiducial model with fragmentation 62 

2.3.4 Influences of the infall model 65 

2.3.5 The radial drift barrier revisited 68 

2.3.6 The fragmentation barrier revisited 74 

11 



Contents 

2.3.7 Dust enhancement inside the snow hne 74 

2.4 Discussion and conclusions 76 

3 The Algorithm 81 

3.1 Numerical schemes 81 

3.1.1 Advection-diffusion Algorithm 81 

3.1.2 Coagulation Algorithm 83 

3.1.3 The fragmentation matrix 85 

3.1.4 Fully implicit scheme for radial motion and coagulation ... 86 

3.2 Test cases 89 

3.2.1 Growth timescales 89 

3.2.2 Radial motion 91 

3.2.3 Temperature structure 91 

4 Dust size distributions in coagulation/fragmentation equilibrium 95 

4.1 Introduction 96 

4.2 Power-law solutions for a coagulation-fragmentation equilibrium . . 97 

4.2.1 The growth cascade 98 

4.2.2 Coagulation fragmentation equilibrium 101 

4.2.3 Fragment dominated regime 102 

4.2.4 Summary of the regimes 103 

4.3 Simulation results interpreted 103 

4.3.1 Constant kernel 104 

4.3.2 Including settling effects 106 

4.3.3 Non-constant kernels 108 

4.4 Grain size distributions in circumstellar disks 109 

4.4.1 Relative velocities 109 

4.4.2 Fragmentation and cratering Ill 

4.4.3 Regime boundaries 113 

4.4.4 Resulting steady-state distributions 115 

4.5 Fitting formula for steady-state distributions 116 

4.5.1 Limitations 116 

4.5.2 Recipe 119 

4.6 Conclusions 121 

4. A Derivation of the fragment dominated size distribution 123 

4.A.1 First integral 124 

4.A.2 Second integral 124 

4.A.3 Third integral 125 

4. A. 4 Deriving the steady-state distribution 126 

5 Dust retention in protoplanetary disks 127 
12 



Contents 



5.1 Introduction 128 

5.2 Model 128 

5.3 Results 131 

5.4 Discussion and conclusions 133 

6 Testing the theory of grain growth and fragmentation by mil- 
limeter observations of protoplanetary disks 135 

6.1 Introduction 136 

6.2 Model description 137 

6.2.1 Disk model 137 

6.2.2 Dust model 138 

6.2.3 Comparison to observations 139 

6.3 Results 140 

6.3.1 Sub-mm fluxes and spectral indices 140 

6.3.2 Radial proflles of the dust opacity index 142 

6.4 Discussion and conclusions 143 

7 Summary and outlook 145 
Bibliography 149 
Acknowledgements 157 



13 



CHAPTER 



1 




Introduction 



The average density in the universe today is about one atom per cubic meter 
(Schneider 2008). Nevertheless, we are hving on an over-density of 29 orders of 
magnitude, known to us as the Earth. Even though the question how our world 
was created is almost as old as mankind itself, we are still far away from being able 
to answer it. Today, astronomers are able to draw a coarse picture of how stars 
and their planetary systems are "produced" but at the same time, many more 
questions have been discovered which have not been posed before. 

This chapter tries to give a very brief overview over the historical aspects and 
the nowadays widely believed scenario of star and planet formation, which sets the 
scene for the following chapters of this thesis. 



1.1 Historical perspective 

The question about the origin of the Earth and our solar system has always been 
particularly attractive as it is one of the keys to constraining the uniqueness of our 
own existence. The earth was believed to be the center of the universe, until 400 
years ago in January 1610, Galileo Galilei pointed his telescope to Jupiter (Drake 
2003). By the discoveries of the moons of Jupiter and the phases of Venus, he 
did not only become the father of modern astronomy but also did he degrade the 
Earth to be just one planet amongst others, orbiting the Sun. This hypothesis was 
soon accepted by other thinkers of his time as evidence was accumulating. 

While previous theories about the world were mostly mystic or philosophical, 
the Age of Enlightenment and the observations of Galilei and Kepler led to more 
mechanical and more realistic philosophies. One of the first of such philosophies 
about the solar system was put forward by Descartes in 1644. He imagined the 
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Figure 1.1: The nebula hypothesis, first suggested by Laplace. In this picture 
of solar system formation, the initial cloud contracts (a) the contraction causes 
the forming disk to spin up due to angular momentum conservation (b). Once 
the centrifugal force equals the gravitational force of the star, the disk fragments 
locally into rings (c) and finally, planets condense out of the rings (d). 



solar system and the universe in total as a system of vortices, thus implying that 
each star itself could host a planetary system like our own. While the idea that the 
solar system was born out of a large vortex was not too far from what we know 
today, his theory still lacked a mathematical basis. He writes (see, Gaukroger 
2006): 

"Finally, we see that, although these whirlpools always attempt a cir- 
cular motion, they practically never describe perfect circles, but some- 
times become too great in width or in length. Thus we can easily 
imagine that all the same things happen to the planets; and this is all 
we need to explain all their remaining phenomena." 

Nevertheless, things were a bit more complicated than Descartes imagined and his 
ideas were quickly swept away by the Newtonian/Keplerian system which seemed 
to explain perfectly the motions of the planets by a mathematical formalism. Yet 
it did not explain the origin of the solar system. 
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1.1 Historical perspective 

During the 18th century, two scenarios were developed, the first of which was 
later known as the tidal theory, first put forward by the french naturalist Georges 
Buffon, the other one was called the nebula hypothesis which was proposed by 
Pierre Simon de Laplace in 1796 (Murdin 2001). In Buffons scenario, a close 
encounter between the sun and a comet should have stripped off material which 
later condenses to form the earth and the other planets. This theory implied 
that the formation of the solar system was a particularly rare event which was 
consistent with the fact that no evidence for other planetary systems was known 
at that time. 

For Laplace, on the other hand, the very uniform motion of the planets and 
their satellites proved that the system was created by a more general law of na- 
ture. Thus he proposed that the solar atmosphere once extended much further out 
through the whole solar system, initially being hot and luminous, and then cooled 
and contracted (see Fig. 1.1). In his view, matter would spin up due to angular mo- 
mentum conservation until it condenses locally, forming rings around the sun out 
of which in turn the planets condense. Other defenders of this scenario include 
Immanuel Kant and Friedrich Wilhelm Herschel who already observed nebulae 
which seemed to resemble the nebulae predicted by this theory. In this scenario, 
planet formation is a natural, universal by-product of star formation which would 
mean that each star could in principle have a planetary system. 

This theory became more and more accepted and few people at that time 
still considered theories of interaction. In 1905, T. C. Chamberlin and F. R. 
Moulton successfully revived Buffons one-century-old dualistic theory (see Brush 
1996), arguing that the Earth could not have formed out of hot gas because its 
atmosphere would have evaporated. They suggested that the solar system formed 
through a close encounter of the sun with another star which induces tidal bulges 
which are then thrown off into the surrounding space. However, unlike in Buffons 
scenario, the planets did not "condense" out of this matter, in their view, but 
were assembled out of small cold particles which they called planetesimals (from 
'infinitesimal planets'). But also the tidal theory had its flaws: Henry N. Russel 
showed that the expelled gas would be far too hot to condense into planets and 
that it is impossible to expel enough material to large enough radii to explain the 
existence of the gas giants (Murdin 2001). 

During the first and the beginning of the second half of the twentieth century, 
many different scenarios and theories were discussed, amongst them the externally- 
triggered cloud collapse model of A. G. W. Cameron, the theory of magnetic 
breaking of Hannes Alfven and Fred Hoyle or the Bondi-Hoyle-Lyttleton cloud 
capture scenario. One of the most influencing ones was by Carl Friedrich von 
Weizsacker who revisited the idea of Descartes' vortices on a more solid basis. He 
was also the first to realize that a turbulent disk would transport mass inside and 
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angular momentum outside (Bodenheimer 2006). 

The subject of star and planet formation experienced an ever increasing im- 
petus since the second half of the last century. Larson (1969), Penston (1969), 
and Shu (1977) derived first solutions of the collapse of self gravitating spheres, 
while the development of computers enabled astronomers to carry out simulations 
of the proposed scenarios which could not be computed analytically. Today, it is 
generally accepted that stars form out of molecular clouds with turbulence being 
both triggering as well as regulating the formation of stars (see, e.g. McKee & 
Ostriker 2007). 

If a cloud core exceeds a certain critical mass, called the Jeans mass^, it be- 
comes gravitationally unstable and will collapse, unless other stabilizing effects, 
such as turbulence, rotational energy or magnetic fields are strong enough to stabi- 
lize the cloud. These effects also gain importance as the cloud collapses: rotational 
energy can cause binary formation and/or the the formation of a circumstellar disk 
(Larson 2003). Additionally, magnetic fields play a role by driving outfiows and 
transporting angular momentum. 

These facts show how diverse and complicated the formation of stars is, al- 
ready on larger scales. Things will not become simpler if we focus on smaller 
scales now. Large parts of our current basic understanding of disk evolution and 
planet formation has to be contributed to a few infiuential papers: Lynden-Bell & 
Pringle (1974) and Shakura & Sunyaev (1973) pioneered the theory of accretion 
disks. These disks are an important aspect of star formation and at the same 
time the birth places of planets. Safronov (1969) and independently Goldreich 
& Ward (1973) envisioned a scenario of planetesimal formation which is (with 
modifications) still discussed today. 

1.2 Theoretical concepts 

Circumstellar disks are the birth places of planets, thus understanding the mech- 
anism which determine their shape and structure is the key to finding out how 
planets can form. The evolution of the gas and the dust in circumstellar disks, 
the growth and fragmentation of the precursors of planet formation and some ob- 
servational implications thereof are the topic of this thesis. In the following, we 
will discuss the most important theoretical concepts of circumstellar disks and 
planet formation and also some observational constrains which together provide 
the framework of this thesis. 



^Thc Jeans mass can very simplistically be derived by equating the thermal and the gravita- 
tional energy of a sphere of constant density. More sophisticated calculations like Jeans' original 
analysis as well as more recent work derive a critical mass with the same dependencies and only 
slightly different coefRcients (e.g., Larson 1985). 
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1.2 Theoretical concepts 

1.2.1 Structure and evolution of circumstellar disks 

Already the most basic observations of our solar system contain valuable hints 
about the formation of stars and planets: most of the mass (about 99.87%) of the 
solar system is contained in the sun, which consists almost entirely of hydrogen and 
helium. In contrast to this, the angular momentum of Jupiter already exceeds the 
angular momentum of the sun by about two orders of magnitude (e.g., Armitage 
2009). The typical angular momentum of a one solar mass molecular cloud core 
(see Goodman et al. 1993) is another 3-4 orders of magnitude larger than the total 
angular momentum of the solar system (assuming a sphere of constant density in 
solid-body rotation). This tells us that which ever mechanism produced our solar 
system out of such a cloud core needs to transport angular momentum outside 
while most mass in transfered inwards. This is known as the angular momentum 
problem of star formation. 

A collapse of a cloud with non- vanishing angular momentum will form a disk. 
This can be understood if one imagines each parcel of gas to orbit around the 
core center maintaining its angular momentum. If the total angular momentum 
of the previously mentioned sphere of constant density is about J ~ 10^*^ g cm^ 
s^^ (assuming a density of n = 10^ cm^^ and the ratio of rotational over gravita- 
tional energy /3 = 0.02), then the parcel of gas with an average specific angular 
momentum / = J/Mq will hit the equator at a radius of about r = ^j^ — 200 AU 
(Ulrich 1976). At the mid-plane, it will collide with a symmetrically moving parcel 
of gas. This will produce a so called accretion shock where the kinetic energy is 
dissipated and thus a disk is formed. In reality, the physics of a cloud core collapse 
is much more complicated because the rotational support can produce multiple 
stars/disks and turbulence and magnetic fields also play an important role dur- 
ing the collapse. The formation of a disk is however a well established outcome, 
also of more realistic simulations of cloud collapse (see, Inutsuka et al. 2010, for 
state-of-the-art simulations). The disks around new-born stars are usually called 
circumstellar disks, accretion disks (implying an inward flux of mass) or protoplan- 
etary disks (implying that planets will form inside of them). Due to the geometry 
of a disk, the vertical sound- crossing time scale is much shorter than the radial 
one, thus it is usually a good approximation to treat the vertical structure to be 
in an equilibrium, while the radial evolution is happening on larger time scales. 

A flrst order vertical structure of the gas disk can be derived by assuming the 
disk to be isothermal (vertical temperature proflle T{z) = const.) and thin (i.e., 
the height above the mid-plane z is small compared to the distance to the star 
r). Then, the vertical component of the stellar gravity needs to be balanced by a 
vertical pressure gradient, 

Cs^ = -Pg^— ^, (1-1) 
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where Cg is the isothermal sound speed, G the constant of gravity, Pg the gas 
density, and M^, the mass of the central star. From this we derive by integration 
a Gaussian density profile 

Pg = Poexpf-^^j , (1.2) 

where we denoted the mid-plane density as po and defined the pressure scale height 
H^ to be the ratio of the sound speed over the Kepler frequency Vt. 

Having a first expression for the vertical structure, we can now focus on the 
radial structure and its evolution. Depending on the desired resolution/complexity, 
this can be done in 1, 2 or 3 dimensions. However, in this thesis, we are interested in 
(parameter studies of) the long-time evolution of these disks and multidimensional 
calculations are computationally too expensive to simulate several million years of 
disk evolution (especially if grain growth is considered as we will discuss later on). 
Therefore - and also to get an analytical feeling of the disk evolution - we will 
now consider the 1+lD approach, where an axisymmetric disk is assumed to be in 
a vertical equilibrium as discussed above. We are thus only interested in the time 
evolution of the surface density, which is the vertically integrated column density, 
defined as 

poo 

Sg(^) = pg{r,z)dz. (1.3) 

J— 00 

If there is viscosity acting within the disk (we will cover this question further 
below), then two neighboring sheets of gas will interact through friction due to 
the differential rotation. The gas is approximately orbiting in Keplerian rotation. 
Therefore, an inner sheet of gas is rotating faster than an outer one and viscous 
friction thus accelerates the outer sheet while the inner sheet is decelerated. This 
mechanism obviously transports angular momentum outside and at the same time 
allows matter to flow inwards. 

The mathematical description of this process is called the theory of accretion 
disks and was first derived by Lynden-Bell & Pringle (1974). They derived the 
radial velocity of the gas due to mass and angular momentum conservation in the 
presence of a viscosity u to be 

t^gas(r) = -^-^^(SgZ/V^). (1.4) 

The evolution of the surface density is then given by the vertically integrated mass 
conservation equation, and can be written as 



5Sg 3 d 
dt r dr 
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1.2 Theoretical concepts 

(for a derivation, see also Pringle 1981). This equation can now be solved numer- 
ically for a given viscosity and initial condition (see Chapter 2 and 3), but it is 
instructive to discuss some general properties of this equation analytically, which 
we will do in the following. 

Substituting x = l-^Jr and u = Sg y^ in Eq. 1.5, we can derive a diffusion 
equation 

du ^(P'u 

with diffusion constant D = Ylvjx^ (assuming a constant v). The diffusion con- 
stant defines the time on which diffusion acts to smooth out gradients in the 
concentration (in this case: Sg). A viscous disk with size r will therefore evolve 
on the viscous timescale given by 

- ^^ ^^ 
f" vise ~pr ~ • V -'- • ' / 

From the observed live time of circumstellar disks, which is about 3 million years 
(see Section 1.3), we can now get constraints on the viscosity of these disks. It 
turns out that the viscous time scale of a 100 AU disk (with typical values for the 
temperature and surface density) is of the order of 10^^ years if we consider only 
molecular viscosity. This is much longer than the age of the universe. To explain 
the observed disk lifetimes, another mechanism of angular momentum transport 
is needed and the anomalous viscosity provided by turbulence is typically believed 
to be solution (although gravitational instability and disk winds may also play a 
role). 

The source of the turbulence is still unknown. Keplerian disks are linearly 
stable to axisymmetric perturbations (e.g., Safronov 1969), but a number of other 
possible sources have been discussed, such as the convective instability (today 
discarded due to the fact that it would transport angular momentum inside instead 
of outside, see Stone & Balbus 1996), the gravitational instability (Toomre 1964; 
Gammie 2001), or the baroclinic instability (Klahr & Bodenheimer 2003; Lesur 
& Papaloizou 2010) to name a few of them. However none of these has been 
as successful as the magnetorotational instability (MRI), discovered by Balbus & 
Hawley (1991), because it draws its energy directly from the Keplerian shear flow 
(Balbus & Hawley 1998). Therefore it has been shown to be a very reliable and 
effective mechanism of angular momentum transport which is effective enough to 
explain the observed accretion rates. Still, it requires two conditions in order to 
operate in circumstellar disks: outwardly decreasing differential rotation and an 
ionization fraction of at least nion/nneutrai ^ 10^^'^ (see Balbus & Hawley 1998). 

While the former condition is always the case in (nearly) Keplerian disks, the 
latter condition may be questionable. A certain level of ionization is provided 
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by cosmic rays, high energy irradiation of the central star, and by decay of ra- 
dioactive nuclei within the disk. However dense enough regions and/or the large 
scale magnetic field line structure may shield some parts of the disk (e.g., Gammie 
1996; Sano et al. 2000; Turner & Drake 2009). In this case, there will be a zone 
where MRI shuts off (called dead zone) surrounded by an MRI- active layer. Since 
the vertically averaged angular momentum transport in layered (i.e., partly-dead) 
regions is not as effective as in active regions, material will accumulate in the dead 
parts of the disk (Gammie 1996). These layered disks can become unstable and 
produce outbursts (thought to explain the FU Orionis objects, e.g. Armitage et al. 
2001) but stationary solutions are also possible (Terquem 2008). 

Coming back to the evolution of the disks surface density evolution, we note 
that there are possible mechanisms providing the necessary anomalous viscosity, 
but it is still unknown which mechanism operates and how effective it really trans- 
ports angular momentum. It has therefore been proven very helpful to parametrize 
our ignorance in the following way, proposed by Shakura & Sunyaev (1973): tur- 
bulent viscosity is provided by eddies with a certain size L and velocity V . The 
largest reasonable values in a circumstellar disk are the scale height H^ and the 
sound speed Cg, respectively. The viscosity can therefore be written as a fraction 
of a product of these quantities, 

u = atCsHp, (1.8) 

where the dimensionless parameter at specifies the amount of turbulence in the 
disk and is therefore usually bound between < at < 1- Typical values of at 
from simulations of MRI-active disks are in the range of 10^^ to a few times 10^^ 
(e.g., Johansen & Klahr 2005; Dzyurkevich et al. 2010). For a given value of at, 
we can now calculate the viscosity if the temperature profile is known. Realistic 
models need to take into account the effects of external irradiation, viscous heat- 
ing, optical depth, and shadowing. Thus multidimensional radiative transfer is 
in principle necessary to derive a realistic temperature structure, however we are 
only interested in some scaling relations of the surface density in this introductory 
chapter and will therefore assume a power-law profile for the temperature 

T(r)ocr"«, (1.9) 

and also for the surface density 

Eg(r)ocr"P. (1.10) 

Postulating a steady-state of the surface density {-gf = 0) and using the power- 
laws for temperature and surface density in Eq. 1.5, we derive 

3 

p + q = --. (1.11) 
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This means that for a given temperature profile, a disk will evolve towards an 
equilibrium power-law with index p. The time scale of this evolution is given by 
Eq. 1.7. From this result, we can also derive the accretion rate of the disk: the 
product rule applied to Eq. 1.4 gives 

^^ ^ 2r Sg ar ^ ' 

Inserting this equation into the definition of the mass accretion rate, 

M(r) = -27rrSgt;g,„ (1.13) 

results in two contributions to the mass flux, 

M(r) =37rSg;/ + 67rr -.^ ■ (1-14) 

We can see now that this mass accretion rate is constant in a steady state disk 
since the product of Eg • z/ is constant with radius and that the mass in a non- 
stationary disk will be rearranged to approach a steady state. Eq. 1.5 can also be 
solved analytically in a less hand-waving way, which can be found in Lynden-Bell 
& Pringle (1974). Their self-similar solution of the disk evolution shows similar 
features. In general, a viscous disk (with power-laws for temperature and viscosity) 
evolves in the following way: the disk profile in the inner parts of the disk (defined 
by a characteristic radius, see Hartmann et al. 1998, for example) follows the same 
power-law derived above, while in the outer parts, the surface density decreases 
exponentially with radius. The mass accretion rate causes the mass of the disk to 
gradually decrease. In order to achieve this, a small fraction of the mass needs to 
"absorb" the angular momentum. The outer parts of the disk therefore spreads 
while the inner parts accrete onto the star. The final outcome of this evolution is 
the state of lowest energy: almost all the mass is in the central object while almost 
all angular momentum is in a small fraction of the mass at a large radius. This 
arrangement is similar to what is observed in our solar system. 

1.2.2 The dynamics of dust 

So far we have only talked about the gas disk evolution. The canonical value of 
the dust-to-gas ratio in the interstellar medium used in the literature is only 1%, 
thus one might be tempted to neglect it. However, not only is this small fraction 
of the total mass the material out of which terrestrial planets, asteroids and the 
cores of giant planets are made, it is also an important ingredient to many aspects 
of circumstellar disks. Almost all the opacity of the disk is provided by grains. 
They do therefore determine most observable properties of disk. Additionally, 
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they are important for chemical surface reaction (e.g., formation of water) or for 
photoelectric heating which strongly influences the temperature in the upper layers 
of the disk (and also the FUV photoevaporation, see Gorti et al. 2009, for example) 
and also for the ionization fraction, which detemins the coupling of the disk to 
the magnetic field and thus its MRI activity (e.g., Sano et al. 2000; Turner & 
Drake 2009). All these facts tell us how important the radial and also the size 
distribution of dust in circumstellar disks is. In this subsection we will introduce 
the most important aspects wich determine the evolution of circumstellar dust. 

As noted above, the dust mass in circumstellar disk is believed to be (initially) 
only about 1% of the gas mass since the disk material stems from the interstellar 
medium. Therefore, the gas is the dynamically dominating phase. This can only 
change if some transport mechanism effectively accumulates the dust relative to the 
gas until the mass fraction approaches unity. We will therefore begin by discussing 
the transport of dust in both the vertical and the radial direction. Analyzing the 
interplay between gas and dust motion is beyond the scope of this chapter (for 
details, see Chiang & Youdin 2010, and references therein). We will rather discuss 
the steady-state solutions of the dust velocity without considering back reactions 
to the gas. 

While the gas reacts to pressure forces, the dust only feels the gravitational 

forces acting on it and the drag force by which it is coupled to the gas motion. 

The drag force is given by 

in 

-Fdrag = Aw, (1.15) 

where m is the mass of the particle, A.u is the particle velocity relative to the gas 
and Ts is the stopping time, which will be discussed in more detail in Chapter 2. 
This coupling to the gas causes dust particles to move towards the region of highest 
pressure. If we consider a particle at a height z above the mid-plane, then its 
vertical motion would be an oscillation around the mid-plane if it were not coupled 
to the gas. The drag force damps this motion and thus the particle will settle 
towards the mid-plane. Equating the vertical component of the gravitational force 
Fg = —mVt\z and the drag force (Eq. 1.15) yields the terminal settling velocity 

Msctt = -TsVlIz. (1.16) 

It turns out to be very practical to define a dimensionless stopping-time of the 
particle called the Stokes number St. It describes the coupling of the particle to 
the gas motion and is in our context defined as the ratio of the stopping time and 
the orbital time, 

St = rs(^k- (1-17) 

This way, particles with different structure or composition but the same St will 
behave dynamically entirely the same. 
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As noted before, the gas disk is most probably turbulent. The distribution of 
dust is thus mixed along with the turbulent gas motion, again depending on the 
coupling to the gas. The vertical distribution of dust is now determined by the 
two counteracting mechanisms of turbulent mixing and vertical settling. We can 
derive a steady-state distribution from the diffusion-advection equation of the dust 
concentration (e.g., Dubrulle et al. 1995) 



dt 



Msctt Pd - -Dd Pg V 



VPg/ 
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Postulating a steady-state [dp^/dt = 0), using Eq. 1.16 and integration results in 
a Gaussian distribution of the dust concentration with a scale height 




hd = Hp\h7- (1-19) 



This leads to the question what the diffusion constant D^ is. Youdin & Lithwick 
(2007) showed that the dust diffusion coefficient is in good approximation given 
by 

Dd = — ^. (1.20) 

l + St" ^ ' 

It is then typically assumed that the gas diffusivity Dg is equal to the viscosity 
z/. Johansen & Klahr (2005) found that in MRI turbulence, radial diffusion is as 
strong as turbulent viscosity, while the vertical diffusivity is slightly weaker. 

The dust-gas coupling also influences the radial distribution of the dust, again, 
by turbulent mixing and "radial settling" towards the star but also by co-moving it 
along with the radial gas flux. These topics will extensively be discussed in Chap- 
ter 2. Especially the radial drift has important implications for planet formation, 
as we will see in the following: larger bodies (i.e., St » 1) are not significantly 
affected by the gas drag (because of their small surface to mass ratio) and thus 
move on Keplerian orbits. In contrast to this, very small particles (St « 1) are so 
tightly coupled to the gas that they are basically frozen-in to the gas. In between 
these two limits, there is a critical regime around a Stokes number of unity: these 
particles are marginally coupled, thus they are not frozen-in to the gas but are 
orbiting at Keplerian velocity. However unlike the much larger bodies, they are 
still significantly infiuenced by the gas drag. The fact that the gas is orbiting at 
sub-Keplerian velocity due to its pressure support induces a relative velocity be- 
tween the dust particles and the gas. The dust particles loose angular momentum 
because of this head wind and consequently spiral inwards. The resulting radial 
velocity can be written as 

St Hp fd\nP\ 

Mdrift - :; ^^72 -w-, Cs, (1.21) 

1 + SV r \dlnr J 
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as derived by Weidenschilling (1977) and Nakagawa et al. (1986) (here, P = p^Cg 
is the gas pressure). For typical disk models, the time scale of this drift can be very 
short, in the order of only a few local orbits (e.g., Brauer et al. 2007). This drift 
poses a fundamental problem for all theories which assume hierarchical growth 
from small dust particles to planetesimals and will be discussed in further detail 
in Chapter 2. The basic ideas of grain growth and planet (-esimal) formation will 
be introduced in the next section. 

1.2.3 From dust to planets 

The growth from sub-yum sized monomers to larger aggregates is inevitably the 
first stage of the formation of planets. Small dust aggregates are sticky because of 
their large surface to mass ratio and inter-molecular binding forces (e.g. van der 
Waals forces or hydrogen bonds). This fact is proven by both experimental (e.g., 
Heim et al. 1999; Poppe et al. 2000) and theoretical work (e.g., Dominik & Tielens 
1997). In order to assemble larger aggregates, these sub-/im sized dust particles 
need to collide. There are several mechanisms which could provide relative motion 
of the dust particles: Brownian motion, the random motion caused thy the thermal 
motion of the gas, leads to random velocities between small particles. The smaller 
the particle, the more it is influenced by Brownian motion. Radial, azimuthal, 
and vertical motion of the dust (as discussed in the previous section) is also size 
dependent, thus, particles with different sizes will have different velocities which 
leads to relative motion between grains of different sizes. Different sized dust 
grains couple to different sized eddies of the turbulent motion of the gas. The 
relative particle velocities induced by turbulence have been derived by Ormel & 
Cuzzi (2007). 

The relative velocities together with the efficient sticking of small dust aggre- 
gates leads to a fast growth of dust; sizes of around centimeters can be reached 
within less than 1000 years at a few AU, as was shown by many authors (Wei- 
denschilling 1980, 1995; Nakagawa et al. 1981; Tanaka et al. 2005; Dullemond & 
Dominik 2005; Brauer et al. 2008a, to name a few). However, two problems arise: 
particles of decimeter sizes at a few AU (or millimeter sizes at around 100 AU) 
have a Stokes number of unity and, once they are formed, are quickly lost to- 
wards the star. Additionally, the relative velocities due to drift and turbulence 
also reach their maximum at these sizes and are large enough to fragment or erode 
the particles (e.g., Brauer et al. 2008a; Blum & Wurm 2008). Additionally, it is 
questionable if the efficient sticking assumption which is often used is still valid 
in this size regime (Youdin 2010). Aggregates of mm size do not seem to stick at 
any velocity, but rather bounce or fragment (Blum & Wurm 2008). This obsta- 
cle might be overcome by the collisions of aggregates of very different sizes (e.g., 
Wurm et al. 2005; Teiser & Wurm 2009). 
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The gravitational instability of the dust is an attractive scenario which basi- 
cally "jumps" over this problematic size regime. In the original idea, developed by 
Safronov (1969) and independently by Goldreich & Ward (1973), growth and verti- 
cal settling of dust particles forms a dust layer at the mid-plane of the disk which 
becomes unstable to its own gravity and thus collapses to form planetesimals. 
Weidenschilling (1980) showed that this scenario does not hold since the thin dust 
layer is unstable to the Kelvin-Helmholtz instability because of the shear at the 
boundary between the dust and the gas. It turns out that mid-plane turbulence is 
both a curse and a cure (Youdin 2010): it may on average stir up concentrations 
of solids, but it also does create random over-densities. The streaming instability 
(Youdin & Goodman 2005) is particularly effective in creating dust over-densities, 
it can even create gravitationally bound clumps, as shown by Johansen et al. 
(2007). However, this scenario still requires particles of around decimeter sizes, 
underlining the importance of the initial "sticking" growth of dust. 

Leaving these problems aside, a new mode of growth is reached, once particles 
have grown to sizes where gravity dominates the accretion of solid bodies (however, 
gas drag can still influence the accretion process, see Ormel & Klahr 2010). The 
gravitational attraction effectively increases the cross section (Lissauer 1993) 






'grav 



\Au) 



a„r.., = Tia" I 1+ f^) ) , (1.22) 



where A-u is the relative velocity of the planetesimals and -Ucsc is the escape velocity 
of the two colliding bodies. The accretion rate in a system where gravitational 
focusing is important (wesc » A-u) scales superlinearly with mass (Ormel et al. 
2010a). Under these conditions, it can be shown (e.g. Armitage 2009) that the ratio 
of the masses of two bodies growth with time. This is called runaway growth and 
will continue until the runaway bodies have reached sizes of a few 100 km (Wetherill 
& Stewart 1989; Thommes et al. 2003; Ormel et al. 2010b). At this point, the 
largest body gravitationally stirs up the distribution of smaller bodies and the 
accretion rate then becomes again sub-linearly with mass (e.g. Ida & Makino 1993; 
Kokubo & Ida 2000; Thommes et al. 2003). Ultimately, the growth of the largest 
body (called the oligarch) is stalled once it has accreted all material within its 
gravitational reach, the so-called feeding zone. The resulting protoplanets typically 
have about the mass of Mars in the terrestrial zone (Thommes & Duncan 2006). 
Unless the disk mass is scaled up by a factor of 10, this scenario fails to explain 
the formation of Earth and Venus and also the formation of the cores of gas 
giant planets in situ: according to Pollack et al. (1996), a core of about 10 M® 
is needed to start runaway gas accretion and thus produce gas giants. Possible 
solutions to this include the migration of small bodies by aerodynamic drag as 
well as the migration of the planets and cores. Scattering between the cores is also 
an attractive scenario as it could possibly explain the formation of Uranus and 
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Neptune which are too massive to be explained by core accretion at their current 
location in the disk (Thommes et al. 2003). 

An entirely different approach to planet formation is the disk instability sce- 
nario: while the disk is typically gravitationally stable in the inner parts, this 
might be debatable in the outer parts where the disk is colder and the orbital fre- 
quency is lower. This can be seen from Toomre's stability criterion (Toomre 1964) 
which states that disk becomes gravitationally unstable if the stability parameter 

Q ^ ^, (1.23) 

becomes smaller than unity. Whether the disk really fragments can however not 
be seen just from this parameter. Before collapse sets in (i.e., 1 < Q < 2), the disk 
develops spiral arms which tend to stabilize the disk by rearranging the surface 
density (Gammie 2001). Additionally, the disk must be able to cool quickly enough 
so that a contracting clump is not sheared apart by the Keplerian rotation. Several 
authors (e.g.. Boss 1997, 2000, 2007; Rice et al. 2005; Meru & Bate 2010) have 
produced very different results which depend highly on the efficiency of cooling in 
the disk (and thus also on the numerical scheme). Today it seems unlikely that disk 
instability alone can explain giant planet formation since the simulations typically 
produce planets at the high end of the mass scale (Rafikov 2005), however some 
of the planets detected by direct imaging (see Sect. 1.3.3) are candidates for this 
formation scenario. 



1.3 Observational constraints 

After having introduced the basic concepts theoretically, it is time to provide 
some proof from the observational side which led to our current understanding 
(and also the problems) of disk evolution and planet formation. In this section, we 
will discuss the most important constraints, derived from observations, ordered by 
size: from observations of gas and tiny dust grains, to the meteorites in our own 
system, and finally to the observed population of planets in and beyond our solar 
system. 

1.3.1 Disk observations 

Circumstellar disks can be detected by several techniques, each tracing specific 
features, temperatures and/or regions of the disk. Combining this knowledge 
enables us not only to proof the existence of disks but also to probe the physics 
which is happening inside of them. 

Most of the disk mass is expected to reside in gas (because disks are made 
out of interstellar material), however this component is the hardest to detect as 
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1000.0 



Figure 1.2: Origin of the spectral energy distribution: at larger wavelength, the 
disk emission outweighs the stellar spectrum. Generally, the shorter wavelength 
emission comes from hotter, inner regions of the disk. Large wavelength (sub- 
millimeter) probe the dust mass in the outer, optically thin regions. Figure adapted 
from Dullemond et al. (2007). 



it is often frozen out, especially in the outer, cold parts of the disk and in the 
mid-plane. 

Another gas feature is a strong Ha excess, which is is believed to be due to an 
accretion flow onto the stellar surface (see Bouvier et al. 2007). Traditionally, T 
Tauri stars are divided into classical T Tauri stars (CTTS) and weak-lined T Tauri 
stars (WTTS), depending on the strength of the emission lines. The accreting 
material has to come from circumstellar material and thus CTTS are thought 
to be younger objects which still have an accretion disk around them. Indeed, 
the accretion rate (typically of the order of 10^^ solar masses per year) and the 
fraction of young stars which have detectable disks (see below) seem to correlate 
(in the same way) with the age of the system (e.g., Calvet et al. 2000; Haisch 
et al. 2001; Sicilia-Aguilar et al. 2006a,b), indicating a typical hfetime of a gaseous 
circumstellar disk to be around 3 million years. Only few gaseous disks have been 
detected around stars older than 6 million years but dust can continue to exist in 
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debris disks, which are probably the remnants of planet (-esimal) formation, similar 
to the asteroid belts in our own solar system. 

The accretion signature and the observed lifetimes of disks show that some 
mechanism must efficiently transport angular momentum outwards and mass in- 
wards. In general, the transition from accreting to non-accreting seems to be very 
rapid, in the order of 10^ years (Simon & Prato 1995; Wolk & Walter 1996; An- 
drews & Williams 2005). This is in contradiction with simple viscous evolution 
which predicts a gradually fading disk (e.g., Hueso & Guillot 2005). Photoevapo- 
rative dissipation of the disk (Hollenbach et al. 1994; Clarke et al. 2001; Alexander 
et al. 2006) or planets blocking the accretion of the disk (e.g., Calvet et al. 2002) 
are some of the frequently suggested scenarios. However at every age some sources 
do not show any signs of accretion. The reason for this is yet unknown. 

Much easier to detect than the gas is the dust which re-radiates the stellar 
irradiation in the infrared (IR). This emission occurs as a long- wavelength excess 
above the stellar photospheric flux, where different wavelength correspond to dif- 
ferent temperatures and thus different regions in the disk. The disk geometry and 
the according spectral energy distribution (SED) is schematically shown in Fig. 1.2. 
In fact this IR excess was the first observational evidence for disks around young 
stars (e.g., Strom et al. 1989; Beckwith et al. 1990). 

The emission of warm dust in the inner disk can be used as a tracer for grain 
growth at smaller sizes: the silicate 10 /im emission feature is prominent for small 
dust sizes (around sub-/im) but disappears as the grain size is increased beyond 
a few /im. van Boekel et al. (2003) found evidence for grain growth at these 
aggregate sizes as one would expect from both theory and laboratory experiments 
(Blum & Wurm 2000; Blum et al. 2000). 

While the near- and mid-IR traces mostly the warm/hot inner regions and the 
surface layers of the optically thick dust disk, observations at millimeter wavelength 
probe the cold outer regions of the disk. Radiation at these wavelength comes from 
the far end of the disk SED, which is close to the Ray leigh- Jeans limit. The disk 
SED can be approximated by (e.g., Beckwith et al. 1990; Andrews & Williams 
2007) 

pF,^PK,B,{T)Mjd\ (1.24) 

if the disk is optically thin and the emission is from an isothermal region. Here 
V denotes the frequency, t^v the opacity (per gram of dust), B^{T) the Planck 
function at a characteristic temperature T, M^ the dust mass, and d the distance 
to the source. Eq. 1.24 shows that under these conditions, the SED is proportional 
to the product of dust mass and opacity. For an assumed opacity law, the dust 
mass can be determined. The total mass of the disk can then roughly be estimated 
of one assumes the dust-to-gas ratio to be the same as in the interstellar case. 
Disk masses were found to be ranging from 0.005 to 0.14 solar masses spread 
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out over sizes of between 20 to several hundred AU (e.g. Andrews et al. 2009). 
The resolution of observations at these wavelength has reached a state where the 
emission can be radially resolved (e.g. Isella et al. 2010). The Atacama Large 
Millimeter Array (ALMA) will be able to resolve the emission down to a few AU 
for the most nearby disks. 

In the framework of viscous disks and a-turbulence (see Eqs. 1.5 and 1.8), 
spatially resolved observations can be used to constrain important parameters of 
the theoretical disk models. Results indicate for example that the turbulence 
parameter should be in the range of at = 5 x 10^^-8 x 10^^ (Andrews et al. 2009), 
which compares well to values predicted by MRI simulations. 

Additionally, the mm-wavelength regime is also an important probe of grain 
growth: Eq. 1.24 shows that the spectral dependency of the SED is given by the 
dust opacity and the Planck function. If we know the mid-plane temperature (or 
are sufficiently deep in the Rayleigh- Jeans regime), measuring the slope of the SED 
(by observations at two or more different wavelengths) provides a direct measure 
of the spectral index of the dust opacity law n^ccX^^ . (3 depends on the dust size 
distribution (and the opacity model). For a power-law dust size distribution, 

n(a)cc|Q else ' ^^'^^^ 

we can now compare which parameters q and Cmax are possible for a certain mea- 
sured value of P (see Fig. 1.3). For example, a /3 value of unity can only be reached 
by a distribution which has a maximum grain size Omax of around a few millimeters. 
Several authors (e.g. Testi et al. 2003; Wilner et al. 2005; Rodmann et al. 2006; 
Ricci et al. 2010) detected emission which indicates mm- or cm-sized dust particles 
at several tens of AU. On the other hand, Brauer et al. (2007) has shown that (for 
reasonable disk parameters) this distribution of grains should drift towards the 
star on a short time scale. A solution to this problem is yet unknown, but parti- 
cle clumping in streaming instabilities or spiral density waves could significantly 
decrease the drift time scale. 

1.3.2 Meteorites 

Additional constrains are coming from a quite different direction, namely mete- 
orites: they provide ample information with a precision unmatched by any astro- 
nomical observation about the formation time scales of solar system bodies. The 
largest class of meteorites, the chondrites are extraordinarily diverse, but they are 
undifferentiated which means that their parent body has not been molten which 
would have differentiated it into different phases. While the oldest known rocks 
on earth are "only" about 4280 million years old (O'neil et al. 2008) (some zircon 
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Figure 1.3: Dust opacity (upper panel) and the resulting opacity index for differ- 
ent dust size distributions. The parameters of the dust size distribution are the 
maximum grain sizes amax and the (negative) power-law index q (see Eq. 1.25). 



minerals are as old as 4404 million years, see Wilde et al. 2001), meteorites can 
be significantly older: the age of calcium- aluminium rich inclusions (CAIs), which 
can be found in chondrites, was determined to be 4568.3 million years, with an 
uncertainty of only 700 000 years (Burkhardt et al. 2008). This way, meteorites are 
the "fossils of the planet formation process" , conserving hints about the structure 
and the chemical composition of the early solar system. 

Most chondrites consist of three main components: Firstly, the CAIs, up to cm- 
sized particles, which are considered to be the first condensates in the solar system. 
Secondly, chondrules, small, mostly spherical silicate crystals with diameters in the 
range of 50 /xm to several mm (Scott 2007) comprise more than 70% by volume 
of the ordinary chondrites (Youdin 2010). Thirdly, the space between CAIs and 
chondrules is filled with the matrix, which consists of sub-/im sized grains of similar 
minerals. A cut through a fragment of the famous Allende meteorite is shown in 
Fig. 1.4. 

The decay of radioactive nuclei can be used to infer both the ages and the 
formation time-scales of these components (for a review, see Russell et al. 2006). 
Results indicate that CAIs are the oldest components and were formed during 
a time interval of about 0.25 million years. Chondrules were formed during an 
extended period of around 2 million years, starting 1-2 million years after the 
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formation of CAIs (Kleine et al. 2009). Today, it is believed that chondrules also 
formed earlier than this, but these early chondrules have been incorporated into 
parent bodies which differentiated due to the heat from ^^Al-decay. Dating of 
iron meteorites indicate that they have formed within < 1 million years after CAI 
formation, they are therefore often thought to be core material of the early parent 
bodies which were differentiated by the ^^Al-heating. 

Despite the ample laboratory data (or perhaps due to it) there is no satisfying 
theory which explains the formation of these bodies. Still two major points that 
have been learned from meteoritic data should be mentioned: firstly, apart from 
some very volatile materials, the composition of chondritic meteorites is strikingly 
similar to the composition of the sun. This tells us that meteorites (and thus the 
planets) formed from the same material as the sun itself. Secondly, the constraints 
on the lifetime of the solar nebula are in very good agreement with the lifetimes 
of circumstellar disks. 




Figure 1.4: A fragment of the AUende meteorite. 



1.3.3 Exoplanets 

The most important observational constraint is that planets do exist, not only 
in our own solar system. The basic concepts behind the most successful planet- 
searching techniques are rather simple: the radial velocity method measures the 
"wobble" of the star induced by the common motion around the center of mass, 
while the transit method uses the dip in the light curve to proof the presence of an 
extrasolar planet. Together, these two methods are responsible for more than 93% 
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of all exoplanet detections. Still, it is the high precision of the measurements that 
is needed, which makes detections of exoplanets such a challenging task. Therefore 
it has taken 40 years after the first proposal of the radial velocity method by Struve 
(1952) that an extrasolar planetary companion (around a neutron star) could be 
detected by Wolszczan & Frail (1992). Detections of planets around solar- type 
stars followed a few years later (Mayor & Queloz 1995). 

Today, the number of extrasolar planets is ever increasing (unlike the number of 
planets in the solar system, which has decreased), proving that planet formation is 
indeed a common by-product of star formation. Extrapolation of recent discoveries 
suggest that about 20% of all stars have planets orbiting within 20 AU (Gumming 
et al. 2008). Population synthesis models predict an even higher rate: about 
30-40% of all FGK stars should host planets (Mordasini et al. 2009). 

Fig. 1.5 shows the mass of the currently ^ known planets as function of their 
semimajor axis. The known population of exoplanets suggests that stars with 
higher metallicity are more likely to host planets (Fischer & Valenti 2005), which 
is in general in agreement with the core accretion scenario. 

Recently, the first direct imaging detections of extrasolar planets have been 
reported (e.g., Kalas et al. 2008; Marois et al. 2008), which have typically several 
Jupiter masses and are located at large distances from the star (hundreds of AU). 
Some of these are thought to have formed by disk instability (Meru & Bate 2010). 



1.4 Summary 

Summarizing all of this, problems with planet formation seem to come in all sizes: 
smaller particles, as discussed above, face the problem that growth can be halted 
by charging of the grains (see Okuzumi 2009), by fragmentation and erosion (see 
Chapter 5), or by bouncing (see Zsom et al. 2010). If larger grains are formed, 
they seem to be lost quickly to the star due to the radial drift. And even if one 
can explain the formation of planetary cores, these bodies still face the problem of 
migration due to interaction with the disk. Additionally, planet formation needs to 
happen relatively quickly (i.e., a few million years), otherwise the gas disk, which 
is needed for the formation of gas giant planets, will have disappeared. Thus, 
we can see that we are able to understand the formation of the solar system and 
the planets far better than Descartes, however we are still far away from a final, 
consistent answer. 



^Data of July 22, 2010 from exoplanet .eu 
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Figure 1.5: The mass of the currently known planets and exoplanets as function 
of their semimajor axis. The solid black circles (•) denote transiting planets, the 
empty circles (o) denote radial velocity detections, and all other detections (pulsar 
timing, microlensing, and direct imaging) are displayed as black dots (•). The red 
symbols represent the solar system planets. 

1.5 About this thesis 

This thesis focuses on the very first stages of planet formation, namely the growth 
of sub-/im particles to planetesimals along with the evolution of the radial distri- 
bution of these particles. Dust is an important ingredient in circumstellar disks, 
not only because terrestrial planets and the cores of giant planets are made out 
of it, but also it is the main source of opacity in disks (thus, important for in- 
terpreting observations), it influences the coupling of the disk to magnetic fields 
by influencing the ionization fraction (important for the viscous evolution of the 
whole disk) and it is an important catalyst in disk chemistry, especially for the 
formation of water. Therefore, we build a numerical code which is able to trace 
the radial and size evolution of dust suspended in a gaseous disk which itself is 
evolving due to viscous transport of angular momentum. 

Chapter 2 describes the physics of growth and radial transport of dust as well 
as the physics of the gas disk which dominates the dynamics of the dust. This is 
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the first work which takes into account the disk build-up phase and subsequent gas 
disk evolution in calculations of the dust evolution. We investigate how the dust- 
to-gas ratio (which is an important parameter for both observations of disks and 
the theory of planetesimal formation) evolves over time if coagulation and frag- 
mentation are taken into account. We also revisit the barriers to planet formation 
by sticking collisions, namely the radial drift and the fragmentation barrier. We 
show that radial variations in the tensile strength of particles (caused by compo- 
sitional changes due to the presence of ices) can cause a strong enhancement of 
dust inside the snow line. The fragmentation velocity turns out to be one of the 
most important parameters determining the particle shape as well as the dust disk 
shape and lifetime. 

The outward spreading of the gas disk is able to carry dust along. Whether and 
to which extent this effect works strongly depends upon the maximum particle size 
and thus also on the critical fragmentation velocity of the gains. There is always a 
radius where even the smallest grains decouple from the spreading gas disk. Future 
observations (e.g., with ALMA) might be able to detect the different outer radii 
of the gas and the dust disks. 

Chapter 3 discusses the numerical implementation of the model. The radial 
transport of the gas and the evolution (in size and radial position) of dust are 
solved using implicit integration schemes. We also show test cases which validate 
the results of our numerical simulations by comparison to analytic estimates. 

Chapter 4 focuses on the grain distribution in the scenario of an equilibrium 
between coagulation and fragmentation. We generalize previous analytical work 
on the slope of steady-state grain size distribution by taking into account both 
fragmentation and coagulation at the same time. These results not only validate 
the distributions we derived in the previous chapters but they also give us the 
opportunity to understand the influence of parameters like the temperature or the 
gas surface density on the size distribution of grains in circumstellar disks. We 
also derive a fltting recipe which can readily be used for further modeling such 
as grain surface chemistry in disks, dust-dependent MRI simulations or radiative 
transfer modeling. 

Chapter 5 investigates the the fact that disks are generally observed to be dusty 
for several million years and how this can be understood if grain fragmentation is 
considered. We show that grain fragmentation at impact velocities similar to what 
is found in experimental studies causes grains to be continuously fragmented to 
small enough sizes to avoid the radial drift regime and to keep the optical depth 
as high as it is observed. 
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Chapter 6 is an apphcation of our previously derived steady state dust size 
distributions to the field of millimeter observations. A grid of different input 
parameters is used to calculate the dust distribution in steady-state. From this, 
we can predict the spectral indices and flux levels at millimeter wavelength and 
compare this to an existing sample of observations. We show that typical values 
can explain the low spectral indices and also the flux levels of a subset of the 
observations. We propose several possibilities which further reduce the fluxes 
while keeping the spectral indices low enough. 

Chapter 7 summarizes our flndings and discusses future prospects of this thesis 
which are planned or already in progress. 
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Evolution of gas and dust in 
protoplanetary disks 

From Birnstiel, Dullemond, & Brauer 2010a, A&A, 513, 79 



Abstract 



Current models of the size- and radial evolution of dust in protoplane- 
tary disks generally oversimplify either the radial evolution of the disk (by 
focussing at one single radius or by using steady state disk models) or they 
assume particle growth to proceed monodispersely or without fragmenta- 
tion. Further studies of protoplanetary disks - such as observations, disk 
chemistry and structure calculations or planet population synthesis models 
- depend on the distribution of dust as a function of grain size and radial 
position in the disk. We attempt to improve upon current models to be 
able to investigate how the initial conditions, the build-up phase, and the 
evolution of the protoplanetary disk influence growth and transport of dust. 
We introduce a new model similar to Brauer et al. (2008a) in which we now 
include the time-dependent viscous evolution of the gas disk, and in which 
more advanced input physics and numerical integration methods are imple- 
mented. We show that grain properties, the gas pressure gradient, and the 
amount of turbulence are much more influencing the evolution of dust than 
the initial conditions or the build-up phase of the protoplanetary disk. We 
quantify which conditions or environments are favorable for growth beyond 
the meter size barrier. High gas surface densities or zonal flows may help 
to overcome the problem of radial drift, however already a small amount of 
turbulence poses a much stronger obstacle for grain growth. 
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Evolution of gas and dust in protoplanetary disks 
2.1 Introduction 

The question of how planets form is one of the key questions in modern astronomy 
today. While it has been studied for centuries, the problem is still far from being 
solved. The agglomeration of small dust particles to larger ones is believed to be at 
least the first stage of planet formation. Both laboratory experiments (Blum et al. 
2000) and observations of the 10 /im spectral region (Bouwman et al. 2001; van 
Boekel et al. 2003) conclude that grain growth must take place in circumstellar 
disks. The growth from sub-micron size particles to many thousand kilometer 
sized planets covers 13 orders of magnitude in spatial scale and over 40 orders 
of magnitude in mass. To assemble a single 1 km diameter planetesimal requires 
the agglomeration of about 10^^ dust particles. These dynamic ranges are so 
phenomenal that one has to resort to special methods to study the growth of 
particles though aggregation in the context of planet (esimal) formation. 

A commonly used method for this purpose makes use of particle size distribu- 
tion functions. The time dependent evolution of these particle size distribution 
functions has been studied by Weidenschilling (1980), Nakagawa et al. (1981), 
Dullemond & Dominik (2005), Brauer et al. (2008a) (hereafter B08a) and others. 
It was concluded that dust growth by coagulation can be very quick initially (in 
the order of thousand years to grow to centimeter sized aggregates at 1 AU in the 
solar nebula), but it stalls around decimeter to meter size due to what is known as 
the "meter size barrier" : a size range within which particles achieve large enough 
velocities to undergo destructive collisions and fast radial inward drift toward the 
central star (Weidenschilling 1977; Nakagawa et al. 1986). 

While the existence of this meter size barrier (ranging in fact from a couple 
of centimeters to tens of meters at 1 AU) has been known for over 30 years, 
a thorough study of this barrier, including all known mechanisms that induce 
motions of dust grains in protoplanetary disks, and at all regions in the disk, for 
various conditions in the disk and for different properties of the dust (such as 
sticking efficiency and critical fragmentation velocity), has been only undertaken 
recently (see B08a). It was concluded that the barrier is indeed a very strong 
limiting factor in the formation of planetesimals from dust, and that special particle 
trapping mechanisms are likely necessary to overcome the barrier. 

However, this work was based on a static, non-evolving gas disk model. It 
is known that over the duration of the planet formation process the disk itself 
also evolves dramatically (Lynden-Bell & Pringle 1974; Hartmann et al. 1998; 
Hueso & Guillot 2005), which may influence the processes of dust coagulation 
and fragmentation. It is the goal of this work to introduce a combined disk- 
evolution and dust-evolution model which also includes additional physics: we 
include relative azimuthal velocities, radial dependence of fragmentation critical 
velocities and the Stokes-drag regime for small Reynolds numbers. 
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The aim is to find out what the effect of dislc formation and evolution is on 
the process of dust growth, how the initial conditions affect the final outcome, 
and whether certain observable signatures of the disk (for instance, its degree of 
dustiness at a given time) can tell us something about the physics of dust growth. 

Moreover, this model will serve as a supporting model for complementary mod- 
eling efforts such as the modeling of radiative transfer in protoplanetary disks 
(which needs information about the dust properties for the opacities) and model- 
ing of the chemistry in disks (which needs information about the total amount of 
dust surface area available for surface chemistry). In this chapter we describe our 
model in quite some detail, and thus provide a basis for future work that will be 
based on this model. 

Furthermore, additional physics, such photoevaporation or layered accretion 
can be easily included, which we aim to do in the near future. 

As outlined above, this model includes already many processes which influence 
the evolution of the dust and the gas disk. However, there are several aspects we 
do not include such as back-reactions by the dust through opacity or collective 
effects Weidenschilling (1997), porosity effects (Ormel et al. 2007), grain charging 
(Okuzumi 2009) or the "bouncing barrier" Zsom et al. (2010). 

This chapter is organized as follows: Section 2.2 will describe all components of 
the model including the radial evolution of gas and dust, as well as the temperature 
and vertical structure of the disk and the physics of grain growth and fragmenta- 
tion. In Section 2.3 we will compare the results of our simulations with previous 
steady-state disk simulations and review the aforementioned growth barriers. As 
an application, we show how different material properties inside and outside the 
snow line can cause a strong enhancement of dust within the snow line. Section 2.4 
summarizes our findings. A detailed description of the numerical method along 
with results for selected test cases can be found in Chapter 3. 

2.2 Model 

The model presented in this work combines a ID viscous gas disk evolution code 
and a dust evolution code, taking effects of radial drift, turbulent mixing, coagu- 
lation and fragmentation of the dust into account. We model the evolution of gas 
and dust in a vertically integrated way. The gas disk is viscously evolving after 
being built up by in-falling material from a collapsing molecular cloud core. 

The radial distribution of grains is subject to gas drag, radial drift, and tur- 
bulent mixing. To which extend each effect contributes, depends on the grain/gas 
coupling of each grain size. By simultaneously modeling about 100-200 different 
grain sizes, we are able to follow the detailed evolution of the dust sub-disk being 
the superposition of all sizes of grain distributions. 
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So far, the evolution of the dust distribution depends on the evolving gas disk 
but not vice versa. A completely self consistently coupled code is a conceptually 
and numerically challenging task which will be the subject of future work. 

2.2.1 Evolution of gas surface density 

Our description of the viscous evolution of the gas disk follows closely the models 
described by Nakamoto & Nakagawa (1994) and Hueso & Guillot (2005). In this 
thesis we shall therefore be relatively brief and put emphasis on differences between 
those models and ours. 

The viscous evolution of the gas disk can be described by the continuity equa- 
tion, 

^j^-lli^,ru,) = S,, (2.1) 

where the gas radial velocity Ug is given by (see Lynden-Bell & Pringle 1974) 

3 d 

u„ = 



Ss-a/t dr 



(Sgz/gv^) . (2.2) 



Eg = '\_^ pg(z)dz is the gas surface density, r the radius along the disk mid-plane 
and z/g the gas viscosity. The source term on the right hand side of Eq. 2.1, denoted 
by Sg can be either infall of material onto the disk or outflow. 

The molecular viscosity of the gas is too small to account for relevant accretion 
onto the star, the time scale of viscous evolution would be in the order of several 
billion years. Observed accretion rates and disk lifetimes can only be explained if 
turbulent viscosity drives the evolution of circumstellar disks. Therefore Shakura 
& Sunyaev (1973) parameterized the unknown viscosity as the product of a velocity 
scale and a length scale. The largest reasonable values for these scales in the disk 
are the pressure scale height Hp 

H, = ^ (2.3) 

and the sound speed Cg. Therefore the viscosity is written as 

z/g = a Cs Hp, (2.4) 

where a is the turbulence parameter and a ^ 1. 

Today it is generally believed that disks transport angular momentum by tur- 
bulence, however the source of this turbulence is still debated. The magneto- 
rotational instability is the most commonly accepted candidate for source of tur- 
bulence (Balbus & Hawley 1991). Values of a are expected to be in the range of 
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10 ^ to some 10 ^ (see Johansen & Klahr 2005; Dzyurkevich et al. 2010). Obser- 
vations confirm this range with higher probabihty for larger values (see Andrews 
et al. 2009). 

If the disk becomes gravitationally unstable, large scale mechanisms of angular 
momentum transport such as through the formation of spiral arms come into 
play. The stability of the disk can be described in terms of the Toomre parameter 
(Toomre 1964) 

Q = -^- (2.5) 

Values below a critical value of Qcr = 2 describe a weakly unstable disk, which 
forms non-axisymmetric instabilities like spiral arms. Q values below unity lead 
to fragmentation of the disk. 

The effect of these non-axisymmetric structures is to transport angular momen- 
tum outward and rearranging the surface density in the disk so as to counteract 
the unstable configuration. This mechanism is therefore to a certain extent self- 
limiting. Values above Qcr are not influenced by instabilities, values below Qcr 
form instabilities which increase Q until the disk is marginally stable again. Our 
model approximates this mechanism by increasing the turbulence parameter a 
according to the recipe of Armitage et al. (2001), 

"«^"» + »»^l|-i|fj;Qjl -i|. (2.6) 

where ao is a free parameter of the model which is taken to be 10^^, unless other- 
wise noted. 

During the time of infall onto the disk, we use a constant, high value of a = 0.1 
to mimic the effective redistribution of surface density during the infall phase which 
also increases the overall stability of the disk. Once the infall stops, we gradually 
decrease the turbulence parameter on a timescale of 10 000 years until it reaches 
its input value. 

Eq. 2.1 is a diffusion equation, which is stiff. This means, one faces the problem 
that the numerical step of an explicit integration scheme goes ccAr^ (where Ar is 
the radial grid step size) which would make the simulation very slow. One possible 
solution to this problem is using the method of implicit integration. This scheme 
keeps the small time scales of diffusion i.e. the fast modes in check. We are not 
interested in these high frequency modes, but they would become unstable if we 
used a large time step. With an implicit integration scheme (see Section 3.1.1) 
the time step can be chosen larger without causing numerical instabilities, thus 
increasing the speed of the computation. 
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2.2.2 Radial evolution of dust 

If the average dust-to-gas ratio in protoplanetary disks is in the order of 10^^ (as 
found in the ISM), then the dust does not dynamically influence the gas, while the 
gas strongly affects the dynamics of the dust. 

Thermally, however, the dust has potentially a massive influence on the gas disk 
evolution through its opacity, but we will not include this in this work. Therefore 
the evolution of the gas disk can, in our approximation, be done without knowledge 
of the dust evolution, while the dust evolution itself does need knowledge of the 
gas evolution. 

There might be regions, where dust accumulates (such as the mid-plane of the 
disk or dead-zones and snow-lines) and its influence becomes significant or even 
dominant but we will not include feedback of such dust enhancements onto the 
disk evolution in this work. 

For now, we want to focus on the equations of motion of dust particles under 
the assumption that gas is the dominant material by mass. The interplay between 
dust and gas can then be described in terms of a dimensionless coupling constant, 
the Stokes number which is defined as 

St = — , (2.7) 

Ted 

where Ted is the eddy turn-over time and Tg is the stopping time. 

The stopping time of a particle is defined as the ratio of the momentum of a 
particle divided by the drag force acting on it. There are four different regimes 
for the drag force which determine the dust-to-gas coupling (see Whipple 1972; 
Weidenschilling 1977). Which regime applies to a certain particle, depends on 
the ratio between mean free path A^fp of the gas molecules and the dust particle 
size a (i.e. the Knudsen number) and also on the particle Reynolds-number Re = 
2aM/z/moi with z/moi being the gas molecular viscosity 

'^mol = i^wAmfp, (2.8) 

u the mean thermal velocity. The mean free path is taken to be 

Amfp = (2.9) 

no-H2 

where n denotes the mid-plane number density and o"h2 = 2 x 10^^^ cm^. 
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The different regimes^ are 



^ for Ajnfp/a > I Epstein regime 



„ ^°" for Re < 1 Stokes regime 1 



QyO.6 1.4 ,,0.4 



(2.10) 



Pgn 



for 1 < Re < 800 Stokes regime 2 
for Re > 800 Stokes regime 3 



Here u denotes tlie velocity of the dust particle with respect to the gas, u = Cs-y/vr/S 
denotes the mean thermal velocity of the gas molecules, ps is the solid density of 
the particles and pg is the local gas density. 

For now, we will focus on the Epstein regime. To calculate the Stokes number, 
we need to know the eddy turn-over time. As noted before, our description of 
viscosity comes from a dimensional analysis. We use a characteristic length scale L^ 
and a characteristic velocity scale Vc of the eddies. This prescription is ambiguous 
in a sense that it does not specify if angular momentum is transported by large, 
slow moving eddies or by small, fast moving eddies, that is 

z/g = (a'^K) • (a^-'^Le). (2.11) 

This is rather irrelevant for the viscous evolution of the gas disk, since all values 
of q lead to the same viscosity, but if we calculate the turn-over-eddy time, we get 

2-kLc ^l-2q 1 {O TO\ 

Ted = —TT- = a " — . (2.12) 

The Stokes number and therefore the dust-to-gas coupling as well as the relative 
particle velocities strongly depend on the eddy turnover time and therefore on q . 
In this work q is taken to be 0.5 (following Cuzzi et al. 2001; Schrapler & Henning 
2004) which leads to a turn-over-eddy time of 

Ted = ^. (2.13) 

The Stokes number then becomes 

St = ^-| fora<^A^fp. (2.14) 



^It should be noted that "Stokes regime" refers to the regime where the drag force on a 
particle is described by the Stokes law - this is not directly related to the Stokes number. 
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The overall radial movement of dust surface density S^ can now be described by 
an advection-difFusion equation, 

aSd 1 d 



+ 



(rFtot)=0, (2.15) 



dt r dr\ 

where the total Flux, Ftot has contributions from a diffusive and an advective flux. 
The diffusive part comes from the fact that the gas is turbulent and the dust 
couples to the gas. The dust is therefore turbulently mixed by the gas. Mixing 
counteracts gradients in concentration, in this case it is the dust-to-gas ratio of 
each size that is being smoothed out by the turbulence. The diffusive flux can 
therefore be written as 

The ratio of gas diffusivity D^ over dust diffusivity Z^d is called the Schmidt num- 
ber. We follow Youdin & Lithwick (2007), who derived 

Sc ^ ^ = 1 + St^ (2.17) 

and assume the gas diffusivity to be equal to the turbulent gas viscosity z/g. 
The second contribution to the dust flux is the advective flux, 

Fadv = Sd-Mr, (2.18) 

where Ur is the radial velocity of the dust. There are two contributions to it, 

1 + St^ ~ St + St" 

The first term is a drag term which comes from the radial movement of the gas 
which moves with a radial velocity of -Ug, given by Eq. 2.2. Since the dust is coupled 
to the gas to a certain extend, the radially moving gas is able to partially drag the 
dust along. 

The second term in Eq. 2.19 is the radial drift velocity with respect to the gas. 
The gas in a Keplerian disk does in fact move sub-keplerian, since it feels the force 
of its own pressure gradient which is usually pointing inwards. Larger dust grains 
do not feel this pressure and move on a keplerian orbit. Therefore, from a point 
of view of a (larger) dust particle, there exists a constant head wind, which causes 
the particle to loose angular momentum and to drift inwards. This depends on 
the coupling between the gas and the particle and is described by the second term 
in Eq. 2.19. «„ denotes the maximum drift velocity of a particle, 

«n = -E^ ■ -^, (2.20) 

z Pg iik 



u, = -^^ ^^^. (2.19) 
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which has been derived by Weidenschilhng (1977). Here, we introduced a radial 
drift efficiency parameter E^. This parameter describes how efficient the radial 
drift actually is, as several mechanisms such as zonal or meridional flows might 
slow down radial drift. This will be investigated in Section 2.3.5. 

Putting all together, the time dependent equation for the surface density of 
one dust species of mass rrii is given by 



dj:i 1 d 
— - -\ 

dt r dr 



^-■"*-^-l(|l-^- 



^g 



SI (2.21) 



where we have included a source term S^ on the right hand side which can be 
positive in the case of infall or re-condensation of grains and negative in the case 
of dust evaporation or outflows. This source term does not include the sources 
caused by coagulation and fragmentation processes (see Section 2.2.6). All sources 
will be combined into one equation later which is implicitly integrated in an un- 
split scheme (see Chapter 3). 

Note that both, the diffusion coefficient and the radial velocity depend on the 
Stokes number and therefore on the size of the particle. 

2.2.3 Temperature and vertical structure 

The vertical structure can be assumed as being in hydrostatic equilibrium at all 
times if the disk is geometrically thin {Hp{r)/r « r) and the vertical sound crossing 
time is much shorter than the radial drift time scale of the gas. The isothermal 
vertical density structure is then given by 

p^{z) = p, exp (^-^ ^y (2.22) 

where 



Po = -^=^. (2.23) 

V 2 TV Hp 

The viscous heating is given by Nakamoto & Nakagawa (1994) 



Q, = S,.,(r^y, (2.24) 

were z/g denotes the turbulent gas viscosity and Q]^ the Kepler frequency. 

Nakamoto & Nakagawa (1994) give a solution for the mid-plane temperature 
with an optically thick and an optical thin contribution, 

Sg a c^ f^k + T;^, (2.25) 
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r [AU] 

Figure 2.1: Total amount of in-fallen surface density as function of radius according 
to the Shu-Ulrich infall model (see Section 2.2.5) assuming a centrifugal radius of 
8 AU (solid line), 33 AU (dashed line), and 100 AU (dash-dotted hue). 



where we used z/g = aCgHp and the approximation tr/p = kr/p Eg/2. «;r and 
Kp are Rosseland and Planck mean opacities which will be discussed in the next 
section. 

Tirr contains contributions due to stellar or external irradiation. Here, we use 
a formula derived by Ruden & Pollack (see Ruden & Pollack 1991, App. B), 



-t irr -^ * 



37r \ r 



^ 2 [ r I \ r 



dlnHp 
dlnr 



(2.26) 



with a fixed dlniJp/dlnr = 9/7, following Hueso & Guillot (2005). 

The main source of opacity is the dust. Due to viscous heating, the temperature 
will increase with surface density. If the temperature rises above ~ 1500 K, the 
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dust (i.e. the source of opacity) will evaporate, which stops the disk from further 
heating until all dust is vaporized. To simulate this behavior in our model, we 
calculate a gas temperature (assuming a small, constant value for gas opacity) 
in the case where the dust temperature rises above the evaporation temperature. 
Then Tmid is approximated by 

-'mid ^'^3'^v-' gas? -^ cvapj) V^'^'j 

only if Tmid from Eq. 2.25 would be larger than Tcvap- 

This is a thermostat effect: if T rises above 1500 K, dust will evaporate, the 
opacity will drop and the temperature stabilizes at T = 1500 K. Once even the 
very small gas opacity is enough to get temperatures > 1500 K, all the dust is 
evaporated and the temperature rises further. 

2.2.4 Opacity 

In the calculation of the mid-plane temperature we use Rosseland and Planck mean 
opacities which are being calculated from a given frequency dependent opacity 
table. The results are stored in a look up table and interpolated during the calcu- 
lations. The opacity table is for a mixture of 50% silicates and 50% carbonaceous 
grains. 

Since these are dust opacities, we convert them from cm? /g dust to crn^/g gas 
by multiplying the values with the dust-to-gas ratio eg, which is a fixed parameter 
in our model, taken to be the canonical value of 0.01. This assumes that the mean 
opacity of the gas is very small and that the dust-to-gas ratio does not change 
with time. To calculate the opacity self-consistently, the total mass of dust and 
the distribution of grain sizes has to be taken into account, meaning that the dust 
evolution has a back reaction on the gas by determining the opacity. For now, our 
model does not take back-reactions from dust to gas evolution into account. Only 
in the case where the temperature rises above 1500 K, the drop of opacity due to 
dust evaporation is considered, as described above. 

2.2.5 Initial infall phase: cloud collapse 

The evolution of the protoplanetary disk also depends on the initial infall phase 
which builds up the disk from the collapse of a cold molecular cloud core. This 
process is still not well understood. First similarity solutions for a collapsing 
sphere were calculated by Larson (1969) and Penston (1969). Shu (1977) found 
a self similar solution for a singular isothermal sphere. The Larson & Penston 
solution predicts much larger infall rates compared to the inside-out collapse of 
Shu {rhui ~ 47 c^/G and 0.975 c^/G respectively). 
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More recent work has shown that the infall rates are not constant over time, 
but develop a peak of high accretion rates and drop to smaller accretion rates 
at later times. The maximum accretion rate is about 13 c^/G if opacity effects 
are included (see Larson 2003). Analytical, pressure-free collapse calculations of 
Myers (2005) show similar behavior but with a smaller peak accretion rate of 
mill = 7.07Cg/G. They also argue that the effects of pressure and magnetic fields 
will further increase the time scales of cloud collapse. 

This initial infall phase is important since it provides the initial condition of 
the disk and also influences the whole simulation by providing a source of small 
grains and gas at larger distances to the star during later times of evolution. 

It should be noted that several groups perform 3D hydrodynamic simulations 
of collapsing cloud cores which show more complicated evolution (e.g., Banerjee 
& Pudritz 2006; Whitehouse & Bate 2006). However, to be able to study general 
trends of the infall phase, we use the Shu-model since it is adjustable by a few 
parameters whose influences are easy to understand. In this model the collapse 
proceeds with an infall rate of mjn = 0.975 c^/G which stays constant throughout 
the collapse. 

We assume the singular isothermal sphere of mass Mdoud to be in solid body 
rotation Vt^. If in-falling material is conserving its angular momentum, all matter 
from a shell of radius r^. will fall onto the star and disk system (with mass ^Ticcnt(^)) 
within the so called centrifugal radius. 



G rncent \t) 



^centr(t) = ^ ^ %,, , (2.28) 



where G is the gravitational constant and r^ = 0.975 • Cg t/2. The path of every 
parcel of gas can then be described by a ballistic orbit until it reaches the equatorial 
plane. The resulting flow onto the disk is 



Sd(r,t) =2pi(r,t)-M,(r,t), (2.29) 



where 



and 



u^{r,t) = \ ju, (2.30) 






as described in Ulrich (1976). Here, fi is given by 

/i = Vl-V^ccntr(r,t). (2.32) 
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The centrifugal radius can therefore be approximated by (cf. Hueso & Guillot 
(2005)) 



W,W ^ 1.4 ( j^)' (^)' (^-^JL_)-\u, p.33) 



We admit that this recipe for the formation of a protoplanetary disk is perhaps 
oversimplified. Firstly, as shown by Visser et al. (2009), the geometrical thickness 
of the disk changes the radial distribution of in-falling matter onto the disk surface, 
because the finite thickness may "capture" an in-falling gas parcel before it could 
reach the midplane. Secondly, star formation is likely to be messier than our simple 
single-star axisymmetric infall model. And even in such a simplified scenario, the 
Shu inside-out collapse model is often criticized as being unrealistic. However, it 
would be far beyond the scope of this chapter to include a better infall model. 
Here we just want to get a feeling for the effect of initial conditions on the dust 
growth, and we leave more detailed modeling to future work. 

2.2.6 Grain growth and fragmentation 

The goal of the model described in this thesis is to trace the evolution of gas and 
dust during the whole lifetime of a protoplanetary disk, including the grain growth, 
radial drift and turbulent mixing. 

Here, the problem arises that radial drift and coagulation "counteract" each 
other in the regime of St = 1 particles: coagulation of smaller sizes restores the 
population around St = 1, whereas radial drift preferentially removes these parti- 
cles. To be able to properly model this behavior, the time step has to be chosen 
small enough if the method of operator splitting is used. 

The upper limit for this time step can be very small. If larger steps are used the 
solution will "flip-flop" back and forth between the two splitted sub-steps of motion 
and coagulation, and the results become unreliable. A method to allow the choice 
of large time steps while preserving the accuracy is to use a non-splitted scheme 
in which the integration is done "implicitly" . B08a already use this technique to 
avoid flip-flopping between coagulation and fragmentation of grains, and we refer 
to that paper for a description of the general method. What is new in the current 
work is that this implicit integration scheme is extended to also include the radial 
motion of the particles. So now radial motion, coagulation and fragmentation are 
done all within a single implicit integration scheme. See Chapter 3 for details. 

Smoluchowski equation 

The dust grain distribution n{m, r, z), which is a function of mass m, distance to 
the star r and height above the mid-plane z, describes the number of particles per 
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Figure 2.2: Sources of relative particle velocities considered in this model (Brow- 
nian motion velocities are not plotted) at a distance of 10 AU from the star. The 
turbulence parameter a in this simulation was 10^^. It should be noted that 
relative azimuthal velocities do not vanish for very large and very small particles. 



cm^ per gram interval in particle mass. This means that the total dust density in 
g cm^^ is given by 



p{r,z 



poo 

)= n[m,r, z) ■ m dm. 

JO 



(2.34) 



With this definition of n(m,r,z), the coagulation/fragmentation at one point 
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in the disk can be described by a general two-body process, 

r) "^ 

—n{m,r,z) = M[m,m',m")x 

° (2.35) 

X n{m', r, z) ■ n{m", r, z) dm dm", 

where M{m,m',m,") is called the kernel. In the case of coagulation and fragmen- 
tation, this is given by 

M{m,, m,', m!') = -K{rri , m!') ■ 5{m! + m!' — m) 



— K{m , m") ■ 5{m" — m) 

-\ — L{m! , m") ■ S{m,, m\ m") 

— L{m' , m") ■ 6{m — m"). 



(2.36) 



For better readability, the dependency of M on radius and height above the mid- 
plane was omitted here. The combined coagulation/fragmentation kernel consists 
of four terms containing the coagulation kernel K, the fragmentation kernel L and 
the distribution of fragments after a collision S. 

The first two terms in Eq. 2.36 correspond to gain (masses m' and m, — m,' 
coagulate) and loss {m, coagulates with m,') due to grain growth. 

The third term describes the fragmentation of masses m and m,', governed 
by the fragmentation kernel L and the fourth term describes the fact that when 
masses m' and m" fragment, they distribute some of their mass via fragments to 
smaller sizes. 

The coagulation and fragmentation kernels will be described in section 2.2.6, 
the distribution of fragments, S, will be described in the next section. 

To be able to trace the size and radial evolution of dust in a combined way, we 
need to express all contributing processes in terms of the same quantity. Hence, 
we will formulate the coagulation/fragmentation equation in a vertically integrated 
way. The vertical integration can be done numerically (as in B08a), however coag- 
ulation processes are most important near the mid-plane, which allows to approx- 
imate the kernels as being vertically constant (using the values at the mid-plane). 
If the vertical dust density distribution of each grain size is taken to be a Gaussian 
(see Section 2.2.6, Eq. 2.51), then the vertical integration can be done analytically, 
as discussed in Section 3.1.2. Unlike the steady-state disk models of B08a which 
have fixed surface density and temperature profiles, we need to recompute the 
coagulation and fragmentation kernels (which are functions of surface density and 
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temperature) at every time step. Therefore this analytical integration also saves 
significant amounts of computational time. 

We therefore define the vertically integrated dust surface density distribution 
per logarithmic bin of grain radius, (T{r, a) as 



a{r,a) = n{r, z,a) ■ m ■ a dz, (2.37) 

J— 00 

where n{r,z,a) and n{r,z,m) are related through m = Air/Spsa^. The total dust 
surface density at any radius is then given by 



Sd(r) = a{r,a) dlna. (2.38) 

Jo 

Defining a{r, a) as in Eq. 2.37 makes it a grid-independent density unlike the 
mass integrated over each numerical bin. This way, all plots of a{r, a) are mean- 
ingful without knowledge of the size grid that was used. Numerically, however we 
use the discretized values as defined in Chapter 3. 

In our description of growth and fragmentation of grains, we always assume 
the dust particles to have a constant volume density meaning that we do not trace 
the evolution of porosity of the particles as this is currently computationally too 
expensive with a statistical code as used in this work. This can be achieved with 
Monte-Carlo methods as in Ormel et al. (2007) or Zsom & Dullemond (2008), how- 
ever these models have do not yet include the radial motion of dust and therefore 
cannot trace the global evolution of the dust disk. 

Distribution of fragments 

The distribution of fragments after a collision, S{m,'m',m"), is commonly de- 
scribed by a power law, 

n{m)dmccm^^dm. (2.39) 

The value C, has been investigated both experimentally and theoretically. Typ- 
ical values have been found in the range between 1 and 2, by both experimental 
(e.g., Blum & Muench 1993; Davis & Ryan 1990) and theoretical studies (Ormel 
et al. 2009). Unless otherwise noted, we will follow B08a by using the value of 

e = 1.83. 

In the case where masses of the colliding particles differ by orders of magnitude, 
a complete fragmentation of both particles is an unrealistic scenario. More likely, 
cratering will occur (Sirono 2004), meaning that the smaller body will excavate a 
certain amount of mass from the larger one. The amount of mass removed from 
the larger one is parameterized in units of the smaller body rUs, 

mont = Xms- (2.40) 
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The mass of the smaller particle plus the mass excavated from the larger one is then 
distributed to masses smaller than rus according to the distribution described by 
Eq. 2.39. In this work, we follow B08a by assuming x to be unity, unless otherwise 
noted. 

Coagulation and fragmentation kernels 

The coagulation kernel K (mi, 7712) can be factorized into three parts, 

K(mi,m2) = Au(mi, 7712) (Tgco{mi, 1712) Pc, (2-41) 

and, similarly, the fragmentation kernel can be written as 

L(mi, 7712) = Au{7ni, 7722) agco(mi, 7722) Pf. (2.42) 

Here, A-u(mi, 7712) denotes the relative velocity of the two particles, o-gco(^^ij ""^2) is 
the geometrical cross section of the collision and pc and pf are the coagulation and 
fragmentation probabilities which sum up to unity. In general, all these factors 
can also depend on other material properties such as porosity, however we always 
assume the dust grains to have a volume density of Ps = 1-6 g cm^^. 

The fragmentation probability is still not well known and subject of both the- 
oretical (Paszun & Dominik 2009; Wada et al. 2008) and experimental research 
(see Blum & Wurm 2008; Giittler et al. 2010). In this work, we adopt the simple 
recipe 

if Au < Uf — 6u 



Pi 



Mf— Am 



if Am > Ui (2.43) 



else 



Su 

with a transition width Su and the fragmentation speed U{ as free parameter which 
is assumed to be 1 m s^^, following experimental work of Blum & Muench (1993) 
and theoretical studies of Leinhardt & Stewart (2009). 

Relative particle velocities 

The different sources of relative velocities considered here are Brownian motion, 
relative radial and azimuthal velocities, turbulent relative velocities and differential 
settling. These contributions will be described in the following, an example of the 
most important contributions is shown in Figure 2.2. 

Brownian motion, the thermal movement of particles, dependents on the mass 
of the particle. Hence, particles of different mass have an average velocity relative 
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to each other which is given by 

8kBT(mi + 1712) ,^ , ,. 

Aubm = \ ^ -■ 2.44 

Y TT mi 1712 

Particle growth due to Brownian relative motion is most effective for small parti- 
cles. 

Radial drift, as described in section 2.2.2 also induces relative velocities since 
particles of different size are differently coupled to the gas. The relative velocity 
is then 

Amrd = \ur{mi) - Ur{m2)\ , (2.45) 

where the radial velocity of the dust, u^ is given by Eq. 2.19. 

Azimuthal relative velocities are induced by gas drag in a similar way as radial 
drift. However while only particles (plus/minus 2 orders of magnitude) around 
St =1 are significantly drifting, relative azimuthal velocities do not vanish for 
encounters between very large and vary small particles (see Figure 2.2). Con- 
sequently, large particles are constantly suffering high velocity impacts of smaller 
ones. According to Weidenschilling (1977) and Nakagawa et al. (1986), the relative 
azimuthal velocities for gas-dominated drag are 



u^ 



Ur, 



1 



1 + stt 1 + st^ 



(2.46) 



where -Un is defined by Eq. 2.20. 

Turbulent motion as source of relative velocities is discussed in detail in Ormel 
& Cuzzi (2007). They also derive closed form expressions for the turbulent relative 
velocities which we use in this work. 

Differential settling is the fifth process we consider contributing to relative 
particle velocities. Dullemond & Dominik (2004) constructed detailed models of 
vertical disk structure describing the depletion of grains in the upper layers of the 
disk. They show that the equilibrium settling speed for particles in the Epstein 
regime is given by -Usett = — -zf^kSt which can be derived by equating the frictional 
force Ffric = —m u/tMc and the vertical component of the gravity force, Fgrav = 
—m z Q"^. To limit the settling speed to velocities smaller than half the vertically 
projected Kepler velocity, we use 

«sett = -z f2k min ( St, - j (2.47) 



for calculating the relative velocities. 
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2.2 Model 



Figure 2.3: Evolution of disk surface density distribution (top) and midplane tem- 
perature (bottom) of the fiducial model described in 2.3.1. 



Since we do not resolve the detailed vertical distribution of particles, we take 
the scale height of each dust size as average height above the mid-plane, which 
gives 

Amds = \hi ■ min(Sti, 1/2) - hj ■ min(Stj, 1/2) | • l^k- (2.48) 

The dust scale height is calculated by equating the time scale for settling, 

z 



t 



sett 



''^sett 



(2.49) 



and the time scale for stirring (Dubrulle et al. 1995; Schrapler & Henning 2004; 
DuUemond & Dominik 2004), 

z^ 

Da 



t 



stir 



(2.50) 



57 



Evolution of gas and dust in protoplanetary disks 

By limiting the vertical settling velocity as in Eq. 2.47 and by constraining the 
dust scale height to be at most equal to the gas scale height, one can derive the 
dust scale height to be 



''' = ''-■ "^" ('' Vmin(St,l/2)(l + St^)j ' ^'-''^ 

This prescription is only accurate for the dust close to the mid-plane, however 
most of the dust (and hence most of the coagulation/fragmentation processes) take 
place near the mid-plane, therefore this approximation is accurate enough for our 
purposes. 



2.3 Results 

2.3.1 Viscous evolution of the gas disk 

We will now focus on the evolution of a disk around a T Tauri like star. We assume 
the rotation rate of the parent cloud core to be 7 x 10^^^ s^^, which corresponds to 
0.06 times the break up rotation rate of the core. The disk is being built-up from 
inside out due to the Shu-Ulrich infall model, with the centrifugal radius being 
8 AU. The parameters of our fiducial model are summarized in Table 2.1. 

Figure 2.3 shows how the gas surface density and the mid-plane temperature 
of this model evolve as the disk gets built up, viscously spreads and accretes 
onto the star. It can be seen that viscous heating leads to a strong increase of 
temperature at small radii. This effect becomes stronger as the disk surface density 
increases during the infall phase. After the infall has ceased, the surface density 
and therefore also the amount of viscous heating falls off. 

This effect also influences the position of the dust evaporation radius, which 
is assumed to be at the radius where the dust temperature exceeds 1500 K. This 
position moves outwards during the infall (because of the stronger viscous heating 
described above). Once the infall stops, the evaporation radius moves back to 
smaller radii as the large surface densities are being accreted onto the star. 

Figure 2.4 shows the evolution of accretion rate onto the star, stellar mass and 
disk mass. The infall phase lasts until about 150 000 years. At this point, the 
disk looses its source of gas and small-grained dust and the disk mass drops off 
quickly until the disk has adjusted itself to the new condition. This also explains 
the sharp drop of the accretion rate. The slight increase in the accretion rate 
afterwards comes from the change in a after the infall stops (see Section 2.2.1). 
Hueso & Guillot (2005) find a steeper, power-law decline of the accretion rate 
after the end of the infall phase because their model does not take the effects of 
gravitational instabilities into account. 
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Figure 2.4: Evolution of disk mass and stellar mass (solid and dashed line on 
left scale respectively) and accretion rate onto the star (dash-dotted line on right 
scale). Adapted from Figure 5 in Hueso & Guillot (2005). 



2.3.2 Fiducial model without fragmentation 

We will now focus on the dust evolution of the disk. This fiducial simulation 
includes only grain growth without fragmentation or erosion. All other parameters 
as well as the evolution of the gas surface density and mid-plane temperature are 
the same as discussed in the previous section. The evolution of this model is 
visualized in Figure 2.5. 

The contour levels in Figure 2.5 show the vertically integrated dust surface den- 
sity distribution per logarithmic bin of grain radius, o"(r, a), as defined in Eq. 2.37. 
The left column of Figure 2.5 shows the results of dust evolution for a steady state 
(i.e. not viscously evolving) gas disk as described in B08a. 

The right column shows the evolution of the dust density distribution of the 
fiducial model, in which the gas disk is gradually built up through infall from 
the parent molecular cloud core, and the gas disk viscously spreads and accretes 
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Figure 2.5: Snapshots of the vertically integrated dust density distributions (de- 
fined in Eq. 2.37) of a steady state disk as in B08a (left column) and of an evolving 
disk (fiducial model, right column). No coagulation is calculated within the evap- 
oration radius (denoted by the dash-dotted line), fragmentation is not taken into 
account in both simulations. The solid line shows the particle size corresponding 
to a Stokes number of unity. Since ast^iccEg (see Eq. 2.14) this curve in fact 
has the same shape as Sg(r), so it refiects, as a "bonus", what the gas disk looks 
like. The radius dividing the evolving disk into accreting and expanding regions 
is marked by the dotted line and the arrows. Particles which are located below 
the dashed line have a positive fiux in the radial direction due to coupling to the 
expanding gas disk and turbulent mixing (particles within the closed contour in 
the upper right plot have an inward pointing fiux). 
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Table 2.1: Parameters of the fiducial model. 



parameter 


symbol 


value 


unit 


turbulence parameter 


a 


10-3 


- 


irradiation angle 


V 


0.05 


- 


cloud mass 


M:loud 


0.5 


Mq 


effect, speed of sound in core 


Cs 


3 X 10^ 


cm s^^ 


rotation rate of cloud core 


^s 


7 X 10-1^ 


s-i 


solid density of dust grains 


Ps 


1.6 


g/cm^ 


stellar radius 


R. 


2.5 


Rq 


stellar temperature 


T. 


4000 


K 



matter onto the star. The solid line marks the grain size corresponding to St =1 
at the given radius. This grain size will be called ast=i hereafter. In the Epstein 
regime, ast=i is proportional to the gas surface density of the disk, which can be 
seen from Eq. 2.14. 

There are several differences to the simulations of grain growth in steady- 
state disks, presented in B08a: viscously evolving disks accrete onto the star by 
transporting mass inwards and angular momentum outwards. This leads to the 
fact that some gas has to be moving to larger radii because it is 'absorbing' the 
angular momentum of the accreting gas. This can be seen in numerical simulations, 
but also in the self similar solutions of Lynden-Bell & Pringle (1974). They show 
that there is a radius R+ which divides the disk between inward and outward 
moving gas. This radius itself is constantly moving outwards, depending on the 
radial profile of the viscosity. 

The radius R+ in the fiducial model was found to move from around 20 AU 
at the end of the infall to 42 AU at 1 Myr and is denoted by the dotted line in 
Figure 2.5. Important here is that small particles are well enough coupled to the 
gas to be transported along with the outward moving gas while larger particles 
decouple from the gas and drift inwards. Those parts of the dust distribution which 
lie below the dotted line in Figure 2.5 have positive fiux in the radial direction due 
to the gas-coupling. 

One might expect that the dotted and dashed lines always coincide for small 
grains as they are well coupled to the gas, however, it can be seen that this is 
not the case. The reason for this is that turbulent mixing also contributes to the 
total flux of dust of each grain size. The smallest particles in between the dotted 
line and the dashed line in the lower two panels of Figure 2.5 do have positive 
velocities, but due to a gradient in concentration of these dust particles, the flux 
is still negative. 

During the early phases of its evolution, a disk which is built up from inside 
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out quickly grows large particles at small radii which are lost to the star by radial 
drift. During further evolution, growth timescales become larger and larger (since 
the dust-to-gas ratio is constantly decreasing) while only the small particles are 
mixed out to large radii. 

The radial dependence of the dust-to-gas ratio after 1 Myr is shown in Fig- 
ure 2.6. These simulations show that the initial conditions of the stationary disk 
models (such as shown in B08a and in the left column of Figure 2.5) are too 
optimistic since they assume a constant dust-to-gas ratio at the start of their sim- 
ulation throughout the disk which is not possible if the disk is being built-up from 
inside out unless the centrifugal radius is very large (in which case the disk would 
probably fragment gravitationally) and grain fragmentation is included to prevent 
the grains from becoming large and start drifting strongly. Removal of larger 
grains and outward transport of small grains lead to the fact that the dust-to-gas 
ratio is reduced by 0.5-1.5 orders of magnitude compared to a stationary model. 
This effect is also discussed in Section 2.3.4. 



2.3.3 Fiducial model with fragmentation 

The situation changes significantly, if we take grain fragmentation into account. 
As discussed in Section 2.2.6, we consider two different kinds of outcomes for grain- 
grain collisions with relative velocities larger than the fragmentation velocity uf. 
cratering (if the masses differ by less than one order of magnitude) and complete 
fragmentation (otherwise). 

For particles larger than about St % 10^^, relative velocities are dominated by 
turbulent motion (and to a lesser extend by vertical settling). Since the relative 
velocities increase with Stokes number (and therefore with grain size), we can 
calculate the approximate position of the fragmentation barrier by equating the 
assumed fragmentation velocity Uf with the approximate relative velocities of the 
particles, 

i^tmax — 7^ o [Z.OZ) 

Particles larger than this size are subject to high-velocity collisions which will 
erode or completely fragment those particles. This is only a rough estimate as the 
relative velocities also depend on the size of the smaller particle and radial drift also 
influence the grain distribution. However Eq. 2.52 reproduces well the upper end of 
the size distribution in most of our simulations and therefore helps understanding 
the influence of various parameters on the outcome of these simulations. 

The evolution of the grain size distribution including fragmentation is depicted 
in Figure 2.7. The initial condition is quickly forgotten since particles grow on very 
short timescales to sizes at which they start to fragment. The resulting fragments 
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Figure 2.6: Comparison of the radial dependence of the dust-to-gas ratio for the 
stationary simulations (thick lines) and the evolving disk simulations (thin lines). 
The four panels compare the results after 1 Myr of evolution with/without frag- 
mentation (left/right column) and with/without effects of non-axisymmetric insta- 
bilities (top/bottom row). It can be seen that the dust-to-gas ratio of the evolving 
disk simulations is almost everywhere lower than the one of the stationary simu- 
lations. See Section 2.3.4 for more details. 



contribute again to the growth process until a semi-equilibrium of growth and 
fragmentation is reached. 

It can be seen that particles stay much smaller than in the model without 
fragmentation. This means that they are less affected by radial drift on the one 
hand and better transported along with the expanding gas disk on the other hand. 
Consequently, considerable amounts of dust can reach radii of the order of 100 AU, 
as seen in Figure 2.8. 

The approximate maximum Stokes number, defined in Eq. 2.52, is inversely pro- 
portional to the temperature (since relative velocities are proportional to Cg), which 
means that in regions with lower temperature, particles can reach larger Stokes 
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Figure 2.7: Evolution of the dust density distribution of the fiducial model as 
in Figure 2.5 but with fragmentation included. The dashed contour line (in the 
lower two panels) around the upper end of the size distribution and around small 
particles at > 60 AU marks those parts of the distribution which have a positive 
(outward pointing) fluxes. 
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numbers. By equating drift and drag velocities of the particles (cf. Eq. 2.19), it 
can be shown that the radial velocities of particles with Stokes numbers larger 
than about q;/2, are being dominated by radial drift. 

Due to the high temperatures below ~5 AU (caused by viscous heating), Stmax 
stays below this value which prevents any significant radial drift within this radius, 
particles inside 5 AU are therefore only transported along with the accreting gas. 
Particles at larger radii and lower temperatures can drift (although only slightly), 
which means that there is a continuous transport of dust from the outer regions 
into the inner regions. Once these particles arrive in the hot region, they get 
"stuck" because their Stokes number drops below a/2. The gas within about 
5 AU is therefore enriched in dust, as seen in Figure 2.8. The dust-to-gas ratio at 
1 AU after 1 Myr is increased by 25% over the value of in-falling matter, which is 
taken to be 0.01. 

Figure 2.8 also shows a relatively sharp decrease in the dust to gas ratio at a 
few hundred AU. At these radii, the gas densities become so small that even the 
smallest dust particles decouple from the gas and start to drift inwards. 

The thick line in Figure 2.8 shows as comparison the dust-to-gas ratio of the 
stationary disk model (cf. left column of Figure 2.5) after 1 Myr, which starts 
with a radially constant initial dust-to-gas ratio of 0.01. 

Figure 2.9 shows the semi-equilibrium grain surface density distribution at 1, 
10 and 100 AU in the fiducial model with fragmentation at 1 Myr. The exact shape 
of these distributions depends very much on the prescription of fragmentation and 
cratering. In general the overall shape of these semi-equilibrium distributions is 
always the same: a power-law or nearly constant distribution (in a) for small 
grains and a peak at some grain size Cmax, beyond which the distribution drops 
dramatically. The peak near the upper end of the distribution is caused by cra- 
tering. This can be understood by looking at the collision velocities: the relative 
velocity of two particles increases with the grain size but it is lower for equal-sized 
collisions than for collisions with particles of very different sizes (see Figure 2.2). 
The largest particles in the distribution have relative velocities with similar sized 
particles which lie just below the fragmentation velocity (otherwise they would 
fragment). This means that the relative velocities with much smaller particles 
(which are too small to fragment the bigger particles but can still damage them 
via cratering) are above this critical velocity. This inhibits the further growth of 
the big particles beyond flmaxj causing a "traffic jam" close to the fragmentation 
barrier. The peak in the distribution represents this traffic jam. 

2.3.4 Influences of the infall model 

In the fiducial model without fragmentation, continuous resupply of material by 
infall is the cause why the disk has much more small grains than compared to 
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Figure 2.8: Evolution of the radial dependence of the dust-to-gas ratio in the fidu- 
cial model including fragmentation with the times corresponding to the snapshots 
shown in Figure 2.7. The initial dust-to-gas ratio is taken to be 0.01. The thick 
dashed curve represents the result at 1 Myr of the static disk model for comparison. 



the stationary disk model (cf. Figure 2.5), which relatively quickly consumes all 
available micrometer sized dust. The effect has already been found in Dominik 
& Dullemond (2008): if all grains start to grow at the same time, then the bulk 
of the mass grows in a relatively thin peak to larger sizes (see Figure 6 in B08a). 
However if the bulk of the mass already resides in particles of larger size, then 
additional supply of small grains by infall is not swept up effectively because of 
the following: firstly, the number density of large particles is small (they may be 
dominating the mass, but not necessarily the number density distribution) and 
secondly, they only reside in a thin mid-plane layer while the scale height of small 
particles equals the gas scale height. 

We studied how much the disk evolution depends on changes in the infall model. 

For a given cloud mass, the so called centrifugal radius rcentr, which was defined 
in Eq. 2.28, depends on the temperature and the angular velocity of the cloud. 
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Figure 2.9: Vertically integrated (cf. Eq. 2.37) grain surface density distributions 
as function of grain radius at a distance of 1 AU (solid), 10 AU (dashed) and 
100 AU (dot-dashed) from the star. These curves represent slices through the 
bottom panel of Figure 2.7. 



Both can be varied resulting in a large range of possible centrifugal radii reaching 
from a few to several hundred AU. Since the centrifugal radius is the relevant 
parameter, we varied only the rotation rate of the cloud core. We performed 
simulations with three different rotation rates which correspond to centrifugal radii 
of about 8 (fiducial model), 33 AU, and 100 AU. For each centrifugal radius, we 
performed two simulations: one which includes effects of gravitational instabilities 
(GI) - i.e. increased a during infall and according to Eq. 2.6 - and one which 
neglects them. 

However for a centrifugal radius of 100 AU, too much matter is loaded onto the 
cold outer parts of the disk and consequently, the disk would fragment through 
gravitational instability. We cannot treat this in our simulations, hence, for the 
case of 100 AU, we show only results which neglect all GI effects. 
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The resulting dust-to-gas ratios are being shown in Figure 2.6. 

Two general aspects change in the case of higher rotation rates: firstly, more of 
the initial cloud mass has to be accreted onto the star by going through the outer 
parts of the disk. Consequently, the disk is more extended and more massive than 
compared to the case of low rotation rate. 

Secondly, as aforementioned the high surface densities in the colder regions at 
larger radii cause the disk to become less gravitationally stable. 

If grain fragmentation is not taken into account in the simulations, both effects 
cause more dust mass to be transported to larger radii. Growth and drift timescales 
are increasing with radius and the dust disk with centrifugal radius of 33 AU 
(100 AU) can stay 5 (35) times more massive than in the low angular momentum 
case after 1 Myr if GI effects are neglected. 

If GI effects are included, matter is even more effectively transported outward, 
the dust-to-disk mass ratio for 8 and 33 AU is increased from 5 to 8. 

However if fragmentation is included, it does not matter so significantly, where 
the dust mass is deposited onto the disk since grains stay so small during the build- 
up phase of the disk (due to grain fragmentation by turbulent velocities) that they 
are well coupled to the gas. Outwards of ~ 10 AU (without GI effects) or of a few 
hundred AU (if GI effects are included), the gas densities become so small that 
even the smallest grains start do decouple from the gas. They are therefore not as 
effectively transported outwards. In these regions, the amount of dust depends on 
the final centrifugal radius while at smaller radii, turbulent mixing quickly evens 
out all differences in the dust-to-gas ratio (see left column of Figure 2.6). 

It can be seen, that in all simulations, the dust-to-gas ratio is lower than in the 
stationary disk model. The trend in the upper right panel in Figure 2.6 suggests 
that for a centrifugal radius of 100 AU and the enhanced radial transport by GI 
effects might have a higher dust-to-gas ratio than the stationary disk model. How- 
ever in this case, the disk would become extremely unstable and would therefore 
fragment. 

The reason for this is the following: to be able to compare both simulations, 
the total mass of the disk-star system is the same as in the stationary disk models. 
How the total mass is distributed onto disk and star depends on the prescription 
of infall. If a centrifugal radius of 100 AU is used, the disk becomes so massive 
that it significantly exceeds the stability criterion Mdisk/^* ^ 0.1. 

2.3.5 The radial drift barrier revisited 

According to the current understanding of planet formation, several mechanisms 
seem to prevent the formation of large bodies via coagulation quite rigorously. The 
most famous ones - radial drift and fragmentation - have already been introduced 
above. Radial drift has first been discussed by Weidenschilling (1977), while the 
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r |AU] 

Figure 2.10: Evolution of the dust surface density distribution of the fiducial model 
at 200 000 years for different drift efficiencies Ed, without fragmentation. The solid 
line denotes the grain size ast=i of particles with a Stokes number of unity. Gas 
outside of the radius denoted by the dotted line as well as particles below the 
dashed line have positive radial velocities. See section 2.3.2 for more discussion. 



importance of the fragmentation barrier (which may prevent grain growth at even 
smaller sizes) was discussed in B08a. In the following, we want to test some 
ideas about how to weaken or to overcome these barriers apart from those already 
studied in B08a. 
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B08a has quantified the radial drift barrier by equating the timescales of growth 
and radial drift which leads to the condition 

^^i(^y*i, (2.53) 

where Cq is the dust-to-gas ratio and t^ and Tg are the drift and growth timescales 
respectively. The parameter 7 describes how much more efficient growth around 
St =1 must be, so that the particles are not removed by radial drift. To overcome 
the drift barrier, obviously either particle growth must be accelerated, or the drift 
efficiency has to be decreased. B08a have numerically measured the parameter 7 
to be around 12. In other words, this means that the growth timescales have to be 
decreased (e.g. by an increased dust-to-gas ratio) until the condition in Eq. 2.53 
is fulfilled. 

However, there are other ways of breaking through the drift barrier. Firstly, 
the drift timescale for St =1 particles also depends on the temperature (via the 
pressure gradient). A simple approximation from Eq. 2.53 with a 0.5 Mq star and 
a dust-to-gas ratio of 0.01 gives 

T<103K {^y\ (2.54) 

which means that particles should be able to break through the drift barrier at 
1 AU if the temperature is below ~100 K (or 200 K for a solar mass star). Dulle- 
mond et al. (2002) have constructed vertical structure models of passively irradi- 
ated circumstellar disks using full frequency- and angle-dependent radiative trans- 
fer. They show that the mid-plane temperature of such a T Tauri like system at 
1 AU can be as low as 60 K. Reducing the temperature by some factor reduces 
the drift time scale by the a factor of similar size which we will call the radial drift 
efficiency ^d (cf. Eq. 2.20). 

Zonal ffows as presented in Johansen et al. (2006) could be an alternative way of 
decreasing the efficiency of radial drift averaged over typical time scales of particle 
growth. Johansen et al. (2006) found a reduction of the radial drift velocity of up 
to 40% for meter-sized particles. 

Meridional ffows (e.g., Urpin 1984; Kley & Lin 1992) might also seem interesting 
in this context, however they do not directly inffuence the radial drift efficiency 
but rather reverse the gas-drag effect. This might be important for small particles 
(which, however are not strongly settling to the mid-plane) but for St =1 particles, 
a would have to be extremely high to have signiffcant inffuence: even a = 0.1 would 
result in a reduction of the particles radial velocity by approximately only a few 
percent. 

Another possibility to weaken the drift barrier is changing its radial depen- 
dence. The reason for this is the following: particle radial drift is only a barrier 
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if it prevents particles to cross the size ast=i- Since particles grow while drifting, 
the particle size corresponding to St =1 needs to increase as well, to be a barrier. 
Otherwise drifting particles would grow (at least partly) through the barrier while 
they are drifting. If ast=i is decreasing in the direction toward the star, then a 
particle that drifts inwards would have an increasing Stokes number even if the 
particle does not grow at all. 



a 10 




Figure 2.11: Dust grain surface density distribution as in Figure 2.10 at 200 000 
years but including the Stokes drag regime. The drift efficiency is set to E^ = 0.75 
and fragmentation is not taken into account. It can be seen that ast=i (solid line) 
increases with radius until about 1 AU, which facilitates the break through the 
drift barrier. 



In the Epstein regime, the size corresponding to St =1 is proportional to the 
gas surface density 



Ost =1 



TTPs 



(2.55) 
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meaning that a relatively flat profile of surface density (or even a profile with 
positive slope) is needed to allow particles to grow through the barrier. However, 
our simulations of the viscous gas disk evolution does not yield surface density 
profiles with positive slopes outside the dust evaporation radius. 

To quantify the arguments above, we have performed simulations with varying 
drift efficiency i?d to test how much the radial drift has to be reduced to allow 
break through. We have additionally included the first Stokes drag regime to see 
how the radial drift of particles is influenced by it. 

Figure 2.10 shows the grain surface density distribution after 200 000 years 
of evolution for three different drag efficiencies. The most obvious changes can 
be seen in the region where the ast^i line (solid line) is relatively fiat: the grain 
distribution is shifted towards larger Stokes numbers. As explained above, particles 
can grow while drifting, which can be seen in the case of E^ = 0.5. The Stokes 
number and size of the largest particles is significantly increasing towards smaller 
radii. However the radial drift efficiency has to be reduced by 80% to produce 
particles which are large enough to escape the drift regime and are therefore not 
lost to the star. 

The situation changes, if the Stokes drag is taken into account: if gas surface 
densities are high enough for the dust particles to get into a different drag regime, 
a change in the radial dependency of ast =i can be achieved. The Epstein drag 
regime for particles with Stokes number of unity is only valid if 

otherwise, drag forces have to be calculated according to the Stokes drag law since 
the Knudsen number becomes smaller than 4/9 (see Weidenschilling 1977). The 
Stokes number is then given by 



v^-KS^n, ^ 1 (2.57) 

9 /imp P ' ^ ' 

with 0"h2 being the cross section of molecular hydrogen. Interestingly, ast=i is 
independent of Sg and proportional to the square root of the pressure scale height 
which decreases towards smaller radii. This means that - as long as the surface 
density is high enough - it does not depend on the radial profile of the surface 
density. In this regime the radial drift itself could move particles over the drift 
barrier since drifting inwards increases the Stokes number of a particle without 
increasing its size. 

Results of simulations which include the Stokes drag are shown in Figure 2.11. 
In this case, particles can already break through the drift barrier if E^ ^ 0.75. This 
value and the position of the breakthrough depends on where the drag law changes 
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Figure 2.12: Dust surface density distributions (top row) and the according solid- 
to-gas ratio (bottom row) for the case of radius-dependent fragmentation velocity 
after 2 x 10^ years of evolution. In the upper row, the vertical dashed line denotes 
the position of the snow line at the mid-plane, the solid line corresponds to ast=i 
and the dot-dashed line shows the approximate location of the fragmentation bar- 
rier according to Eq. 2.52. The snow line on the right plot lies further outside 
since viscous heating is stronger if a is larger. In the bottom row, the solid line 
denotes the vertically integrated dust-to-gas ratio while the dashed line denotes 
the dust-to-gas ratio at the disk mid-plane. The icy dust grains outside the snow 
line are assumed to fragment at a critical velocity of 10 m s^^ while particles inside 
the snow line fragment at 1 m s^^. The plots on the left and right hand side differ 
in the amount of turbulence in the disk [a = 10^^ and 10^^, respectively). 



from Epstein to Stokes regime and therefore on the disk surface density. As noted 
above, larger surface densities generally shift the position of regime change towards 
larger radii making it easier for particles to break through the drift barrier. 

It should be noted that the physical way to avoid the Stokes drag regime is 
to decrease the surface densities, however we chose the same initial conditions for 
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both cases - with and without Stokes drag - and just neglected the Stokes drag in 
the latter computations to be able to compare the efficiency factors independent 
of other parameters such as disk mass or temperature. 

2.3.6 The fragmentation barrier revisited 

In the previous section, we have shown that several mechanisms could allow par- 
ticles to break through the radial drift barrier, however fragmentation puts even 
stronger constraints on the formation of planetesimals. 

As shown by Ormel & Cuzzi (2007), the largest relative velocities are of the 
order of 

Awmax - \/2q;Cs. (2.58) 

If particles should be able to break through the fragmentation barrier, then they 
need to survive these large relative velocities, meaning that A-Umax has to be smaller 
than the fragmentation velocity of the particles, or 



> 



'2a. (2.59) 



This condition is hard to fulfill with reasonable fragmentation velocities, unless 
a is very small. E.g., for a = 10^^ and a temperature of 200 to 250 K, the 
fragmentation velocity needs to be higher than 4 m s^^, which could already seen 
in the simulations by Brauer et al. (2008b), who have simulated particle growth 
near the snow-line. 

Even in the case of very low turbulence, relative azimuthal velocities of large 
(St > 1) and small grains (St « 1) are of the order of 30 m s^^, which means that 
large particles are constantly being 'sand-blasted' by small grains. The only way 
of reducing these velocities significantly is decreasing the pressure gradient (see 
Equations 2.20 and 2.46). 

Another possibility to overcome this problem would be if high-velocity impacts 
of smaller particles would cause net growth, as has been found experimentally by 
Wurm et al. (2005) and Teiser & Wurm (2009). 

Taken together, these facts make environments as the inner edge of dead zones 
ideal places for grain growth (see Brauer et al. 2008b; Kretke & Lin 2007): shutting 
of MRI leads to low values of a, which are needed to reduce turbulent relative 
velocities and the low pressure gradients prevent radial drift and high azimuthal 
relative velocities. 

2.3.7 Dust enhancement inside the snow Une 

As we will discuss in Chapter 5, significant loss of dust by radial drift can be 
prevented by assuring that particles stay small enough and are therefore not in- 
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fluenced by radial drift. For typical values of a (10^^ — 10^^), this means that the 
fragmentation velocity must be smaller than about 0.5-5 m s^^. If particles have 
higher tensile strength, they can grow to larger sizes which are again affected by 
radial drift. 

Typical fragmentation velocities for silicate grains determined both theoreti- 
cally and experimentally are of the order of a few m s^^ (for a review, see Blum 
& Wurm 2008). The composition of particles outside the snow-line is expected 
to change due to the presence of ices. This can influence material properties and 
increase the fragmentation velocity (see Schafer et al. 2007; Wada et al. 2009). 

We have performed simulations with a radially varying fragmentation speed. 
We assume the fragmentation speed inside the snow-line to be 1 m s^^, outside the 
snow-line to be 10 m s^^. It should be noted that we do not follow the abundance 
of water in the disk or the composition of grains, we only assume particles outside 
the snow line to have larger tensile strength due to the presence of ice. To be able 
to compare both simulations, we used the same 1 Myr old 0.09 Mq gas disk around 
a solar mass star as initial condition. The gas surface density profile of this disk 
was derived by a separate run of the disk evolution code. We used this gas surface 
density profile and a radially constant dust-to-gas ratio as initial condition for 
the simulations presented in this subsection. Apart from the level of turbulence, 
the input for both simulations is identical, the results are therefore completely 
independent of uncertainties caused by the choice of the infall model. 

Results of the simulations are shown in Figure 2.12. A one order of magnitude 
higher fragmentation velocity causes the maximum grain size to be about two 
orders of magnitude larger, which follows from Eq. 2.52 since StmaxQcamax in the 
Epstein regime (all particles in these simulations are small enough to be in the 
Epstein regime). This effect can be seen in Figure 2.12. Particles outside the 
snow-line are therefore more strongly drifting inwards (because they reach larger 
Stokes numbers) where they are being pulverized as soon as they enter the region 
within the snow-line. 

Strong drift outside the snow line and weaker radial drift inside the snow line 
cause the dust-to-gas ratio within the snow line to increase significantly (see bottom 
row of Figure 2.12): in the case of a = 10^^, the dust-to-gas ratio reaches values 
between 0.39 and 0.10 in the region from 0.2 to 1.9 AU. 

Simulations for a less massive star (0.5 Mq) show the same behavior, however 
the increase in dust-to-gas ratio is not as high as for a solar mass star (dust-to-gas 
ratio of 0.27-0.20 from 0.2-4 AU). 

This effect is not as significant in the case of stronger turbulence, where the 
maximum dust-to-gas ratio is around 0.027. The evolution of the dust-to-gas ratio 
at a distance of 1 AU from the star is plotted for both cases in Figure 2.13 (the 
minor bump is an artifact of the initial condition). 
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The reason for this difference hes in the locations of the drift and fragmentation 
barriers. The approximate position of the fragmentation barrier is represented by 
the dot-dashed hne in Figure 2.12. The radial drift barrier cannot be defined 
as sharply, however radial drift is strongest at a Stokes number of unity, which 
corresponds to the solid line in Figure 2.12. An increase of a by one order of 
magnitude lowers the fragmentation barrier by about one order of magnitude in 
grain size. 

In the lower turbulence case, the fragmentation barrier lies close to ast=i- Most 
particles are therefore drifting inwards before they are large enough to experience 
the fragmentation barrier. Hence, the outer parts of the disk are significantly 
depleted in small grains. 

In contrast to this case, fragmentation is the stronger barrier for grain growth 
throughout the disk in the high turbulence simulation (right column in Figure 2.12). 
It can be seen that particles of smaller sizes are constantly being replenished by 
fragmentation. 

With these results in mind, the evolution of the disk mass (bottom panel of 
Figure 2.14) seems counter-intuitive: the mass of the high turbulence dust disk is 
decaying faster than in the low turbulence case. This can be understood by looking 
at the global dust-to-gas ratio of the disks (top panel of Figure 2.14) which does 
not differ much in both cases. This means that the increased dust mass loss in the 
high turbulence disk is due to the underlying evolution of the gas disk. Particles 
in the high turbulence disk may have smaller Stokes numbers (causing drift to be 
less efficient), however the inward dragging by the accreting gas is stronger in this 
case. 

To show how much the dust evolution depends on the fragmentation velocity, 
we included the case of a lower fragmentation velocity throughout the disk in 
Figure 2.14. It can be seen that the dust mass is retained at its initial value for 
much longer timescales. 



2.4 Discussion and conclusions 

We constructed a new model for growth and fragmentation of dust in circum- 
stellar disks. We combined the (size and radial) evolution of dust of B08a with 
a viscous gas evolution code which takes into account the spreading and accre- 
tion, irradiation and viscous heating of the gas disk. The dust model includes 
the growth/fragmentation, radial drift/drag and radial mixing of the dust. We 
re-implemented and substantially improved the numerical treatment of the Smolu- 
chowski equation of B08a to solve for the combined size and radial evolution of 
dust in a fully implicit, un-split scheme (see Chapter 3). In addition to that, we 
also included more physics such as relative azimuthal velocities, radial dependence 
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Figure 2.13: Dust-to-gas ratio at a distance of 1 AU from the central star as a 
function of time for the case of low [a = 10^^, solid line) and high [a = 10^^, 
dashed line) turbulence. 



of fragmentation critical velocities and the Stokes drag regime. The code has been 
tested extensively and was found to be very accurate and mass-conserving (see 
Section 3.2). 

We compared our results of grain growth in evolving protoplanetary disks to 
those of steady state disk simulations, similar to B08a. In spite of many differences 
in details, we confirm the most general result of B08a: radial drift and particle 
fragmentation set strong barriers to particle growth. If fragmentation is included 
in the calculations, then it poses the strongest obstacle for the formation of plan- 
etesimals. Very low turbulence [a < 10^^) and fragmentation velocities of more 
than a few m s^^ are needed to be able to overcome the fragmentation barrier in 
the case of turbulent relative velocities. 

This model includes also the initial build-up phase of the disk, which is still a 
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Figure 2.14: Evolution of the global dust-to-gas ratio (top panel) and the dust 
disk mass (bottom panel) for the simulations shown in Figure 2.12. The solid and 
dashed lines correspond to the low and high turbulence case, respectively. The 
dash-dotted line shows the evolution of the disk if a low fragmentation velocity is 
assumed throughout the disk. 
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very poorly understood phase of disk evolution. We use the Shu-Ulrich infall model 
which represents a strong simplification. However, the following novel findings of 
this work do not or only weakly depend on the build-up phase of the disk: 

• Apart from an increased dust-to-gas ratio (as in B08a), other mechanisms 
such as streaming instabilities or a decreased temperature may be able to 
weaken this barrier by decreasing the efficiency of radial drift. We found 
that in simulations without fragmentation the radial drift efficiency needs to 
be reduced by 80% to produce particles which crossed the meter-size barrier 
and are large enough to resist radial drift. 

• If the gas surface density is above a certain threshold (in our simulations 
about 140 g cm^^ at 1 AU or 330 g cm^^ at 5 AU, see Eq. 2.56) then the 
drag force which acts on the dust particles has to be calculated according 
to the Stokes drag law, instead of the Epstein drag law. The drift barrier 
in this drag regime is shifting to smaller sizes for smaller radii (independent 
on the radial profile of the surface density) which means that pure radial 
drift can already transport dust grains over the drift barrier or at least to 
larger Stokes numbers even without simultaneous grain growth. In this case, 
the drift efficiency needs to be reduced only by about 25% to produce large 
bodies. 

• If relative azimuthal velocities are included, then grains with St > 1 are 
constantly 'sand-blasted' by small grains (if they are present) which (in our 
prescription of fragmentation) causes erosion and stops grain-growth even 
in the case of low turbulence. Only decreasing the radial pressure gradient 
significantly weakens both relative azimuthal and radial velocities. Low tur- 
bulence and a small radial pressure gradient together are needed to allow 
larger bodies to form. These conditions may be fulfilled at the inner edge of 
dead zones (Brauer et al. 2008b; Kretke & Lin 2007; see also Dzyurkevich 
et al. 2010). Future work needs to investigate the disk evolution and grain 
growth of disks with dead zones. Our prescription of fragmentation and ero- 
sion may also need rethinking since Wurm et al. (2005) and Teiser & Wurm 
(2009) find net-growth by high velocity impacts of small particles onto larger 
bodies. 

• Higher tensile strengths of particles outside the snow-line allows particles 
to grow to larger sizes, which are more strongly affected by radial drift. 
Particles therefore drift from outside the snow-line to smaller radii where they 
fragment and almost stop drifting (since their radial velocity is decreased by 
almost two orders of magnitude). This can cause an increase the dust-to-gas 
ratio inside the snow-line by more than 1.5 orders of magnitude. 
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• The critical fragmentation velocity and its radial dependence (and to a lesser 
extent the level of turbulence) is a very important parameter determining if 
the dust disk is drift or fragmentation dominated. A drift dominated disk 
is significantly depleted in small grains and only a fragmentation dominated 
disk can retain a significant amount of dust for millions of years as is observed 
in T Tauri disks. 

The following results depend on the build-up phase of the disk. However unless 
the collapse of the parent cloud is not inside-out or so fast to cause disk fragmen- 
tation, we expect only slight alteration of the results: 



• 



• 



• 



Disk spreading causes small particles (< 10 /im) to be transported outward 
at radii of ~60-190 AU even in 1 Myr old disks. 

Small particles provided by infalling material are not effectively swept up by 
large grains if the bulk of the dust mass has already grown to larger sizes. 

In an inside-out build-up of circumstellar disks, grains growth is very fast 
(time scales of some 100 years) because densities are high and orbital time 
scales are small. Large grains are quickly lost due to drift towards the star 
if fragmentation is neglected. Fragmentation is firstly needed to keep grains 
small enough to be able to transport a significant amount of dust to large 
radii by disk spreading and secondly to retain it in the disk by preventing 
strong radial inward drift. 



80 



CHAPTER 
The Algorithm 

Based on Birnstiel, Dullemond, & Brauer 2010a, A&A, 513, 79, Appendix 

3.1 Numerical schemes 

In the next two sections, we will first discuss how the equations of radial evolution 
of gas (Eq. 2.1) and dust (Eq. 2.21) as well as the coagulation/fragmentation 
of dust (Eq. 2.35) are solved separately. In Section 3.1.4 we will then describe 
how this model treats the radial and the size evolution of dust in an un-split, fully 
implicit way. 

3.1.1 Advection-diffusion Algorithm 

To be able to also model both, the evolution of dust and gas implicitly, we con- 
structed a scheme which solves a general form of an advection-diffusion equation. 




dM d , , , 


d 


\ d 


a+a.<-^-"'- 


dx 


h-D^- — 

ox 



-^) 



K + L-M (3.1) 



which can be adapted to both, Eq. 2.1 and Eq. 2.21 by proper choice of parameters 
M, Dd, 5', h, K and L. 

We use a fiux-conserving donor-cell scheme which is implicit in M. The time 
derivative in Eq. 3.1, written in a discretized way becomes 

dN A/'*+^ - A/"* 

— = 1^^^^ ^ (^3^2) 

St ti+i — ti 
where i denotes time-dimension and n denotes space-dimension. 
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The advective part is discretized as 
where F^~^\ and 5'„,i are the future flux and the surface between cell n and cell 

n+5 "-+2 

n + 1 and V„ is the volume of cell n. The advective interface fluxes can then be 
written as 

^It' = A4-max(0,M„+i)+A/;+i-min(0,M„+i) (3.4) 

""•"2 2 2 

F*+\ = Ar„_i • max(0,M„_i) + Afn ■ min(0,M„_i) (3.5) 

"2 2 2 

The diffusive interface flux becomes 



^Mi = ^dn+1 V^ ^^""^^ ^"^ (3-6) 

rf,n+5 a, ".+ 2 "+2 X 1 — T 



CU = ^d,n-i /^n-i '" ! "l'^ (3.7) 



3n A/" gn-1 \/- 

I "-n— 1 

2 X^ X|2— 1 

Working out these equations and separating the values of Af leads to 

K = An- A/^+l + 5n • Ar;+1 + a, ■ Ui,X\ + D„ (3. 

with the coefficients 

/ D, l-h I -S 1 -On-l 

"■ V„ V V ' "-2' "^2 (a;„-a;„_i)/i„_i 

Bn = I- AtLn + ^ i max(0, m„+i ) • 5„+i 
-min(0,M„_i)-5„_i 

D, , 1 -h , 1 -S ,1 -On D , I -h I -S 1 -On 



(3.9) 






£>» = -t\t-K„. 



'■n- 



Eq. 3.8 can now be solved by any matrix-solver, but since it is a tri-diagonal 
matrix, the fastest analytical way to solve it is by forward elimination/backward 
substitution. 

It should be noted that Eq. 3.1 is implicit only in M which means that the 
equations we solve are only implicit in the surface density. In the case of the 
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viscous accretion disk, described by Eq. 2.1, we face the problem that also the 
turbulent gas viscosity z/g depends on the temperature which (in the case of viscous 
heating) depends on the surface density. This can cause numerical instabilities to 
develop. 

To stabilize the code, we use a scheme which estimates the temperature in 
several predictor steps. The actual time step is then done with the predicted 
temperature. 

3.1.2 Coagulation Algorithm 

Discretizing Eq. 2.35 on a mass grid rrii gives 

d 

-~-nk{r,z) = ^Mijkni{r, z) nj{r, z), (3.10) 

where the dust particle number density is 

ni{r,z)= n{'m,r,z)dm (3-11) 

If we assume that the coagulation and fragmentation kernels are constant in z and 
that the vertical distribution of grains is a Gaussian with a scale height according 
to Eq. 2.51, 

nk{r,z)= /^,^. . -exp (- /^ ) , (3.12) 

we can vertically integrate the coagulation/fragmentation equation. 

d f^ 

-A^fc(r) = Ukdz 

J— 00 



(3.13) 



where 



dt 

= y Miik Ni(r) NJr) ^—— x 

= 5]M,,,iV,(r)iV,.(r), 

ij 



Mijk = , Mijk- (3.14) 

l27Jh^^) 
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Discretizing Eq. 3.13 also in radial direction and rewriting it in terms of the 
quantity 

MM = Nj,^-i)-ri (3.15) 



yields 



|a4z = 2^ -M, -AT,, := Su. (3.16) 



where the vector S = {Ski}k=i,m is the source function for each of the m mass bins. 
The numerical change of Mki within a time step At = ti — tj_i is AMm = 
■Sflf^ — Mil- The time-discretized version of Eq. 3.16 then becomes (omitting 
second order terms) 

^ = S ^ (M. + AMO (a;, + aa;,) 

'' J (3.17) 

Since the first term on the right hand side is the explicit source function, we can 
write 



AAfki ^ , ^Mijk + Mjik 

— r— - Jkl + 2_i 

At -^ ri 



S,,+Y,}^bl]l±^^.M,,.ANa. (3.18) 



*j 



Using the vectors AN = {AJ\f]k^i,nm ^^^ ^ = {A/'}fc=i,n„, this can be rewritten 
in matrix notation, 

■ 1 



Were 



^^ J)AN = S. (3.19) 



J.. = 1] ^'^' ;" ^''' Af,i (3.20) 

denotes the Jacobian of the source function and 1 the unity matrix. The solution 
for the future values can now be derived by inverting the matrix in Eq. 3.19, 

^i+i = N^ + AN 

• / 1 \"' (3.21) 
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3.1.3 The fragmentation matrix 

The discretization of the coagulation/fragmentation equation contains some pit- 
falls, which are mostly related to numerical rounding errors: already the growth 
from 0.1 /im to just a decimeter sized object spans a mass range of 18 orders of 
magnitude, which is beyond the numerical precision of a double precision floating 
point number. Problems arise for example if a large body grows by accumulating 
very small grains. The numerical scheme should thus be written in a way where 
addition of very different sized numbers or subtraction of almost-equal sized num- 
bers are avoided. In Brauer et al. (2008a), the coagulation matrix was rewritten 
in order to conserve the total mass, in this subsection, we will show, how the 
fragmentation matrix needs to be handled. 

The fragmentation part of the Smoluchowski equation is written in a de- 
scretized way as 
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(3.22) 



where Lij is the fragmentation kernel, Skij is the distribution of fragments to bin 
A;, which are produced by collisions of rrii and nij. The fraction {rrii + mj]/mk 
converts the source terms from number densities to masses since we want to con- 
serve the total mass. We again face the problem of numerically subtracting large, 
almost equal numbers. To be as accurate as possible, it is best to subtract these 
numbers analytically, if possible, since the computer variables we use have "only" 
16 decimals of precision. To be able to do that, we have to rewrite Eq. 3.22 keeping 
in mind that the fragmentation matrix is symmetric {Lij = Lji) 
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(3.23) 



Now in our prescription of cratering, masses w-big and mbrcak collide, w-break excavates 
X times its own mass from rribig and fragments itself if it is one order of magnitude smaller 
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in mass. The total mass of fragments is mfrag = (1 + x) ^'T'brcak; which is distributed 
according to the power-law described in Eq. 2.39 while the larger mass looses x^^brcak) 
which leaves 

W-rcst = "Ibig - X^T-break- (3.24) 

Numerically, mrest lies somewhere between rribig and the next smaller bin, which we will 
denote as mbig-i- The grid distance is defined as 

Am = TTtbig — 'TZbig-i. (3.25) 

Therefore, rTii-cst will be redistributed into the adjacent grid points according to the 
Podolak prescription. We need to take care since if x or the grid resolution is is large, 
then rrircst might be smaller than rribig-i- For typical cases of x and the grid resolution, 
this is rarely the case. (1 — e) • mrcst will be added to rribig and e • rrircst will be added to 
"T-big-i) where e is defined as 

e = ^^^1B££^. (3.26) 

Am ^ ^ 

Now in the cases where m^ = mbig we can write the contributions to m-big, normalized 

to mtotal = "T-big + W-brcak aS 

'TT-rcst '^big 



Skij = (1 — e)- 

''"■total "Hotal 
"vmi i_ — f^m ^ — mu:„ 

(3.27) 



"T-total "T-total 
"T-big - X"ibreak " e"lrcst " "^big 



'T'ltotal 
'^brcak '^rest 

= -X e • 

"T-totla "T-total 

By separating of this special case in the code, the mass is conserved several orders of 
magnitude better than without rewriting it. 

3.1.4 Fully implicit scheme for radial motion and 
coagulation 

To be able to solve the radial motion and the Smoluchowski equation fully implicitly, 
we rewrite Equation 3.8 as 

M • N'^+^ = N^ - D, (3.28) 

where M is the tri-diagonal matrix with entries A, B and C and source term D which 
are given by Eq. 3.9. M is now rewritten by separating off the diagonal terms and by 
dividing by At, 

M = At • M + 1 (3.29) 

which brings Eq. 3.28 in a form similar to Eq. 3.19, 

1 



+ M ) • AN = -M • N* - D/At (3.30) 
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Figure 3.1: Pictographic representation of the matrix on the left hand side of 
Eq. 3.32 with six radial and five mass grid points. White elements are zero, dark 
grey elements contain contributions from radial transport of dust (J), light gray 
elements contain contributions from the coagulation/fragmentation (M), and black 
elements contain both contributions. The upper left and lower 5 rows represent the 
boundary conditions where coagulation/fragmentation is not taken into account. 
The matrix in the simulations would typically have a size of 15000x15000. 



The coagulation/fragmentation is now determined by the matrix J and the source vector 
S and similarly, the radial motion is determined by matrix M and source vector 



R = -M • N* - D/At. 
This allows us to combine both operators into one matrix equation, 



1 
At 



+ M-J-AN = R + S. 



(3.31) 



(3.32) 



In principle, to solve for the vector AN, the matrix on the left hand side of Eq. 3.32 has 
to be inverted. This is a numerically challenging task since the inverse matrix in our 
simulations would have about 150-500 million elements. The matrix on the left hand 
side of Eq. 3.32, however is a sparse matrix (schematically depicted in Figure 3.1). 
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Figure 3.2: Comparison between simulation result and analytical solution for 
growth of equal sized particles. The solid lines denote the position of the peak 
of the grain size distribution. In the top panel, only Brownian motion is con- 
sidered as source of relative particle velocities, in the bottom panel, turbulent 
relative velocities are considered as well. The parameters of these simulation are 



T = 196 K, ps=1.6 g/cm , Sg = 18 g/cm , Sj = 0.18 g/cm and 



a 



10 



-3 



We can therefore solve Eq. 3.32 by an iterative incomplete LU decomposition solver 
for sparse matrices provided by the Sixpack library^ of S. E. Norris. 



^available from www . engineers . aucklcoid . ac . nz/~snor007 
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3.2 Test cases 

3.2.1 Growth timescales 

To check if the numerical implementation presented above accurately solves Eq. 3.32, 
we compare results of the simulation to some analytical solutions: The growth rates 
of particles can be approximated if we assume the grain size distribution to be a delta 
function and the sticking probability to be unity. The increase of mass per collision is 
then given by the mass of the particle, m divided by the time between two collisions, r, 
which can be written as 

^ = "^. (3.33) 

at T 

Using dm = 47r/9sa^da and r = m/{pA o Au), we derive 

^ = ^ An, (3.34) 

where Am is the relative velocity, p^ the dust density, p^ the solid density of the dust 
particles and a = Aira^ the cross section of the collision. 

With this formula, we can estimate the growth time scales if the relative velocities 
of the dust grains is given. For equal sized particles and Brownian motion, we get 

/16 kkT 
Ambm = A/ , (3.35) 

for turbulent motion, the relative velocities are (see Ormel & Cuzzi 2007) 

V2aSt for St « 1 

for St » 1 ■ P-^*"* 




^-(t-to)+aoM (3-37) 



exp(§|^-(t-to)) forSt<l 
^ • (t - to) + ao for St > 1 



ao ■ 

2a 



(3.38) 



for relative velocities due to turbulent motion. 

A comparison between analytical solution and simulation result for Brownian motion 
growth is shown in the top panel of Figure 3.2. It can be seen that the position of the 
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Figure 3.3: Test case for only radial drift, without coagulation. The solid line 
denotes the mass-averaged position of the radial distribution of grains for each 
grain size after 10^ years. The dashed line is the expected solution for a single 
particle, the dotted line denotes the initial radius of the particles. 



peak of the grain size distribution (solid line in Figure 3.2) follows the analytical result 
ofEq. 3.37. 

A similar comparison for the case of relative velocities due to turbulent motion is 
not as straight-forward since both the turbulent relative velocities and the dust scale 
height (Eq. 2.51) are sub-divided into several regimes. We therefore integrated Eq. 
3.34 numerically for the case of Brownian motion and turbulent relative velocities; the 
results are shown in the bottom panel of Figure 3.2. As before, we see that - after the 
initial condition is overcome - the simulation result follows closely the mono-disperse 
approximation. For grains larger than St = 1, the analytical solution is also plotted in 
the bottom panel of Fig. 3.2, almost coinciding with the simulation results. 
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3.2.2 Radial motion 

The radial motion of dust particle was tested in a similar fashion: starting from a grain 
distribution at a radius of 5.5 AU, we let the particles drift (taking the gas drag and the 
radial drift into account, see Eq. 2.19) without coagulation. We compare the results to 
results of a numerical integration of the equation of motion for a single particle. The 
results are shown in Fig. 3.3. 

We find that the size distribution behaves as expected: small particles are well 
coupled to the gas, they almost retain their initial position since the radial motion due 
to gas drag is in the order of 0.01 AU. Larger particles, having a larger Stokes number 
drift towards the star on shorter timescales. Particles of a few centimeters (corresponding 
to St =1) are already lost after about 700 years. 

The mass in all test cases was found to be conserved on the order of 10^^^% of the 
initial value. 

3.2.3 Temperature structure 

We also compared the results of our temperature equation (cf. Eq. 2.25) to a non- 
isothermal vertical disk profile, which will be described in the following. In the non- 
isothermal case, the source terms are z-dependent: 

vise/ \ "^ "'b-' K^) / \ o 

q+ [z) = -a p{z) 17k 

C"" W = -p(.) „P W j^ . exp (-^) 

Where r is the Rosseland optical depth and if the flaring angle, which is usually also 
dependent on radius. Here we take c/? to be a constant with the value 0.05. The factor 
of 1/2 in Eq. 3.39 comes from the fact that one half of the radiation is radiated away 
from - the other half into the disk (see Dullemond et al. 2002). 

Given the source terms, we can solve for the intensity by using the moment equations 
and the Eddington approximation: if the intensity Ii, is nearly isotropic, 

/(r, /i) = a(r) + 6(r) // (3.40) 

the moment equations become 



2J-1 



J = Idfi = a (3.41) 

2 j-i 

1 f^ b 

K = - \ I- fi^dfi = ^ (3.43) 

2 J-i 



3 



which leads to the closure relation 



K=-J (3.44) 
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Figure 3.4: Vertical Structure at 2 AU with a surface density of 94 g cm^^. Shown 
is the temperature profile (top left panel), the heat sources (top right panel), the 
gas density structure (bottom left) and the Rosseland optical depth (bottom right). 



and from the second moment follows 

dK 

Since (see DuUemond et al. 2002, A. 9) 



H 



IdJ 



(3.45) 






47r ' 



(3.46) 



H can be calculated by integrating over the sources starting from the boundary condition 
HiQ) = 0. Knowing H, we can integrate Eq. 3.45 to get J. The necessary boundary 
condition at r = can be determined in the two-stream-approximation (see Rybicki 
&: Lightman 1991). The temperature follows from Dullemond et al. (2002), Eq. A. 13, 
where we use kp{{J7:/(t-q)^'^) instead of kj. It should also be noted that the last term 
in A. 13 of Dullemond et al. (2002) should be ^/(47rp). 
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Figure 3.5: Difference between the approximate mid-plane temperature, calcu- 
lated from Eq. 2.25 and the non-isothermal vertical structure, as described in this 
subsection. 



With a given temperature, the vertical hydrostatic equilibrium can be calculated. 
From P = pc^ follows 

d{pT) _ _ 2 ^"T-P 



dz 



-pz^k 



kh 



(3.47) 



which can be done numerically by using pi = 1 and by evaluating the right hand side as 
the average value between the grid centers. The resulting profile can then be normalized 
to a given surface density. 

The above procedure has to be done iteratively since the temperature profile depends 
on the density profile and vice versa. We repeat the T and p calculation until the 
temperature at any point does not change by more than 1%. An exemplary result for 
the vertical structure of the disk is shown in Fig. 3. 4. An comparison between the results 
of Eq. 2.25 and this non-isothermal case is shown in Fig. 3.5. 
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CHAPTER 




Dust size distributions in 

coagulation/fragmentation 

equilibrium 

FroTn Birnstiel, Ormel, & DuUemond 2010, submitted to A&A 



Abstract 

Grains in circumstellar disks are believed to grow by mutual collisions 
and subsequent sticking due to surface forces. Many fields of research involv- 
ing circumstellar disks, such as radiative transfer calculations, disk chem- 
istry, magneto-hydrodynamic simulations and others largely depend on the 
resulting grain size distribution. 

As detailed calculations of grain growth and fragmentation are both nu- 
merically challenging and computationally expensive, we aim to find simple 
recipes and analytical solutions for the grain size distribution in circumstel- 
lar disks for a scenario in which grain growth is limited by fragmentation 
and radial drift can be neglected. 

We generalize previous analytical work on self-similar steady-state grain 
distributions. Numerical simulations are carried out to identify under which 
conditions the grain size distributions can be understood in terms of a com- 
bination of power-law distributions. A physically motivated fitting formula 
for grain size distributions is derived using our analytical predictions and 
numerical simulations. We find good agreement between analytical results 
and numerical solutions of the Smoluchowski equation for simple shapes of 
the kernel function. The results for more complicated and realistic cases can 
be fitted with a physically motivated "black box" recipe presented in this 
chapter. 
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4. Dust size distributions in coagulation/fragmentation equilibrium 
4.1 Introduction 

Dust size distributions are a fundamental ingredient for many astrophysical models in the 
context of circumstellar disks and planet formation: whenever dust is present, it domi- 
nates the opacity of the disk, thereby influencing the temperature and consequently also 
the vertical structure of the disk (e.g., Dullemond & Dominik 2004). Small grains effec- 
tively sweep up electrons and therefore strongly affect the chemistry and the ionization 
fraction (also via grain surface reactions, see Vasyunin et al., in prep.) and thereby also 
the angular momentum transfer of the disk (e.g., Wardle &: Ng 1999; Sano et al. 2000). 

Today, it is well established that the dust distributions in asteroid belts and de- 
bris disks are governed by a so-called "collision cascade" (see Williams &: Wetherill 
1994): larger bodies in a gas free environment exhibit high velocity collisions (> km 
s^^), far beyond their critical fragmentation threshold which lead to cratering or even 
complete shattering of these objects. The resulting fragments in turn suffer the same 
fate, thus producing ever smaller grains down to sizes below about a few micrometers 
where Poynting- Robertson drag removes the dust particles (e.g., Wyatt et al. 1999). 
The grain number density distribution in the case of such a fragmentation cascade has 
been derived by Dohnanyi (1969) and Williams &: Wetherill (1994) and was found to 
follow a power-law number density distribution n(m) octtt,^" with index a = ^ (which is 
equivalent to n(a)oca^^'^), with very weak dependence on the mechanical parameters of 
the fragmentation process. Tanaka et al. (1996), Makino et al. (1998) and Kobayashi &: 
Tanaka (2010) showed that this result is exactly independent of the adopted model and 
that the resulting slope a is only determined by the mass-dependence of the collisional 
cross-section if the model of collisional outcome is self-similar. The value of -g- agrees 
well with the size distributions of asteroids (see Dohnanyi 1969) and of grains in the 
interstellar medium (MRN distribution, see Mathis et al. 1977; Pollack et al. 1985) and 
is thus widely applied, even at the gas-rich stage of circumstellar disks. 

However, in protoplanetary disks gas drag damps the motions of particles. Very 
small particles are tied to the gas and, as a result, have a relative velocity low enough to 
make sticking feasible. The size distribution therefore deviates from the MRN power-law. 
Theoretical models of grain growth indicate that particles can grow to sizes much larger 
than a few /im (see Nakagawa et al. 1981; Weidenschilling 1980, 1984, 1997; Dullemond 
& Dominik 2005; Tanaka et al. 2005; Brauer et al. 2008a; Birnstiel et al. 2009). Indeed, 
the observational evidence suggests that growth up to cm-sizes is possible (e.g., Testi 
et al. 2003; Natta et al. 2004; Rodmann et al. 2006; Ricci et al. 2010). 

However, as particles grow, they become more loosely coupled to the gas. This results 
in an increase in the relative velocity between the particles, a common feature of most 
sources of particles' relative velocity (i.e., turbulence, radial drift, and settling motions). 
For this reason, we expect that the perfect sticking assumption will break down at some 
point and other collisional outcomes, e.g., bouncing, erosion, and catastrophic disruption 
become possible (see Blum & Wurm 2008; Giittler et al. 2010). It is expected, then, that 
at a certain point growth will cease for the largest particles in the distribution. Collisions 
involving these particles result in fragmentation, thus replenishing the small grains. 
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On the other hand mutual collisions among small particles still result in coagulation. 
As a result, a steady-state emerges. This situation is referred to as a fragmentation- 
coagulation equilibrium. 

The situation in protoplanetary disks differs, therefore, from that in debris disks. In 
the latter only fragmentation operates. The mass distribution still proceeds towards a 
steady state but, ultimately, mass is removed from the system due to, radiation pressure 
or Poynting-Robertson drag. 

In this chapter, we consider the situation that the total mass budget in the system is 
conserved. For simplicity, we will ignore motions due to radial drift in this study. This 
mechanism effectively removes dust particles from the disk as well as providing particles 
with a large relative motion. However, the derived presence of mm-size dust particles in 
protoplanetary disks is somewhat at odds with the usual (laminar) prescriptions for the 
radial drift rate (Weidenschilling 1977). Recently, it was shown that mm observations 
of protoplanetary disks can be explained by steady-state size distributions if radial drift 
is inefficient (Birnstiel et al. 2010b). On the other hand, if radial drift would operate 
as effectively as the laminar theory predicts, then the observed populations of mm-sized 
particles at large radii cannot be sustained (Brauer et al. 2007). Possible solutions to 
reduce the drift rate include bumps in the radial pressure profile (see Brauer et al. 2008b; 
Cossins et al. 2009) or zonal flows (see Johansen et al. 2009). 

In this work we analytically derive steady-state distributions of grains in the presence 
of both coagulation and fragmentation. The analytical predictions are compared to 
numerical simulations and applied to grain size distributions in turbulent circumstellar 
disks. Both theory and simulation results presented in this work are used to derive a 
fitting formula for steady-state grain size distributions in circumstellar disks. 

This chapter is outlined as follows: in Sect. 4.2, we briefly summarize and then gen- 
eralize previous results by Tanaka et al. (1996) and Makino et al. (1998). In Sect. 4.3, we 
test our theoretical predictions and their limitations by a state-of-the-art grain evolution 
code (see Birnstiel et al. 2010a and Chapter 3). Grain size distributions in circumstellar 
disks are discussed in Sect. 4.4. In Sect. 4.5, we present a fitting recipe for these dis- 
tributions that can easily be used in models where grain properties are important. Our 
findings are summarized in Sect. 4.6. 



4.2 Power- law solutions for a 

coagulation- fragmentation equilibrium 

In this section, we begin by summarizing some of the previous work on analytical self- 
similar grain size distributions on which our subsequent analysis is based. We will then 
extend this to include both coagulation and fragmentation processes in a common frame- 
work. Under the assumption that the relevant quantities, i.e., the collisional probability 
between particles, the distribution of fragments, and the size distribution, behave like 
power-laws, we will solve for the size distribution in coagulation-fragmentation equi- 
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librium. For simplicity, we consider a single monomer size only of mass niQ. This is 
therefore the smallest mass in the distribution. 

Another key assumption of our analytical model is that we assume the existence of 
a sharp threshold mass nif, above which collisions always result in fragmentation and 
below which collision always result in coagulation. As explained above, the physical 
motivation for this choice is that relative velocities increase with mass. This also means 
that in our model we neglect collision outcomes like bouncing or erosion. 

In our analysis, we will identify three different regimes, which are symbolically illus- 
trated in Fig. 4.1: 

• Regime A represents a case where grains grow sequentially (i.e. hierarchically) 
by collisions with similar sized grains until they reach an upper size limit and 
fragment back to the smallest sizes. The emerging power-law slope of the size 
distribution depends only on the shape of the collisional kernel. 

• Regime B is similar to regime A, however in this case the fragmented mass is redis- 
tributed over a range of sizes and, thus, influencing the out-coming distribution. 

• In Regime C, the upper end of the distribution dominates grain growth at all 
sizes. Smaller particles are swept up by the upper end of the distribution and 
are replenished mostly by redistributed fragments of the largest particles. The 
resulting distribution function depends strongly on how the fragmented mass is 
distributed after a disruptive collision. 

For each of these regimes, we derive the parameter ranges in which they apply and also 
the slopes of the resulting grain size distribution. 

4.2.1 The growth cascade 

The fundamental quantity that governs the time-evolution of the dust size distribution 
is the collision kernel Cmi,m2- It is defined such that 

Cmi,m2 ■ n(mi) ■ n{m2) dmidm2 (4.1) 

gives the number of collisions per unit time per unit volume between particles in mass 
interval [mi, mi + drrii] and [m,2,m,2 + dm2], where n(m) is the number density dis- 
tribution. Once specified, it determines the collisional evolution of the system. In the 
case that the number density n(m) does not depend on position, Cmi,m2 is simply the 
product of the collision cross section and the relative velocity of two particles with the 
masses mi and m,2. 

We use the same Ansatz as Tanaka et al. (1996), assuming that the collision kernel 
is given in the self-similar form 

Cmum2 =ni\h f ^) , (4.2) 

\mij 
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Case B) 



Case C) 



log(mass) 



sequential growth 



''^o^\^^fragmentation to monomers/^% 



sequential growth 



fragnient distribution 
sweep-up growth 



fragment distribution 





Figure 4.1: Illustration of the three different regimes for which analytical solutions 
have been derived. Case A represents the growth cascade discussed in Sect. 4.2.1, 
case B the intermediate regime (see Sect. 4.2.2) and case C the fragmentation 
dominated regime which is discussed in Sect. 4.2.3. Particles are shattered once 
they reach the fragmentation barrier rrif since collision velocities for particles > rrif 
exceed the fragmentation threshold velocity. 
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where h is any function which depends only on the masses through the ratio of 1712/1711 
and ly is called the index of the kernel or the degree of homogeneity. As we will see in 
the following, v is one of the most important parameters determining the resulting size 
distribution. Different physical environments are represented by different values of i'. 
Examples include the constant kernel (i.e., i^ = 0, mass independent), the geometrical 
kernel (i.e., i' = 2/3, velocity independent) or the linear kernel (i.e., i^ = 1, as for grains 
suspended in turbulent gas). 

It is further assumed that the number density distribution of dust particles follows 
a power- law, 

n{m)=A-m-°'. (4.3) 

The time evolution of the mass distribution obeys the equation of mass conservation, 

dmn(m) dF(m) „ ,, ,, 

^ + -^^ = 0' (^•^) 

ot cm 

where the flux F{m) does not represent a flux in a typical continuous way since coagula- 
tion is non-local in mass space (each mass can interact with each other mass) but rather 
an integration of all growth processes which produce a particle with mass greater than 
m out of a particle that was smaller than m (i.e. collisions of mi < in with any other 
mass m2 such that rrii + m2 > ra). The flux in the case of pure coagulation was derived 
by Tanaka et al. (1996) 

dmi dm2miCmi,m2'n{7ni)n{m2), (4.5) 

mo Jm—mi 

where we changed the lower bound of the integration over nii to start from the mass 
of monomers niQ (instead of 0) and the upper bound of the integral over ?n,2 to a finite 
upper end rrif (instead of infinity) . 

Substituting the definitions of above and using the dimensionless variables xi = 
rrii/m, X2 = m2/m, xq = rrio/m, and Xf = m{/m one obtains 

F{m) = m'^-^a+s r ^^^ r Ax2xl+^-'^x^'^h (—\ , (4.6) 

^ V ' 

:= K 

where K integrates to a constant. 

Postulation of a steady state (i.e. setting the time derivative in Eq. 4.4 to zero), 
leads to the condition 

F(m)ocm''"2a+3 ^ const., (4.7) 

from which it follows, that the slope of the distribution is 

"-^. (4.8) 

This result was already derived for the case of fragmentation by Tanaka et al. (1996) 
and Dohnanyi (1969) and for the coagulation by Klett (1975) and Camacho (2001). 
The physical interpretation of this is a "reversed" fragmentation cascade: instead of a 
resupply of large particles which produces ever smaller bodies, this represents a constant 
resupply of monomers which produce ever larger grains (cf. case A in Fig. 4.1). 
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4.2.2 Coagulation fragmentation equilibrium 

As Tanaka et al. (1996) and Makino et al. (1998) pointed out, the previous result is inde- 
pendent of the model of collisional outcomes as long as this model is self-similar (Eq. 4.2). 
However this is no longer the case if we consider both coagulation and fragmentation 
processes happening at the same time. 

We will now consider the case with a constant resupply of matter, due to the particles 
that fragment above mf. We assume these fragments obey a power-law mass distribution 
and are produced at a rate 

nf(m) = N -mr^, (4.9) 

where ^ reflects the shape of the fragment distribution. 

If there is a constant flux of particles F{m) as defined above, then the flux of frag- 
menting particles (i.e., the flux produced by particles that are growing over the fragmen- 
tation threshold) is given by 

F(mf) = i^ • mf^-2a+3 (4_;^Q) 

where K is the integral defined in Eq. 4.6. 

The resulting (downward) flux of fragments Fi{m) can then be derived by inserting 
Eq. 4.9 into the equation of mass conservation (Eq. 4.4), 

— £!__:. = _rn . rif. (4.11) 

dm 

Integration from monomer size rriQ to m yields 

Ff(m) = -N ■ -^— ■ (n?-^ - rn^^\ . (4.12) 

The normalization factor N can be determined from the equilibrium condition that the 
net flux vanishes, 

Ff(mf) = -F(mf), (4.13) 



and was found to be 



2-5 2-5 

In Eq. 4.12, we can distinguish two cases: 



N={2-i)K^^^ 5-. (4.14) 



rrij 



• If ^ > 2, the contribution of ra^ dominates the term in brackets. This means that 
most of the fragment mass is redistributed to monomer sizes and the situation is 
the same as in the pure coagulation case (cf. Case A in Fig. 4.1). The steady-state 
condition F{m) + F{{m) = (i.e., the net flux is zero) yields that 

is constant, which leads to the same condition as Eq. 4.8. Intuitively, this is clear 
since the majority of the redistributed mass ends up at m ~ ttiq. 



101 



4. Dust size distributions in coagulation/fragmentation equilibrium 

• If .^ < 2, the m-dependence dominates the term in brackets in Eq. 4.12 and 
postulation of a steady-state, 

K-m''-'^'^+^ ^N-^--m'^-^, (4.16) 



leads to an exponent 



--•^^^. ("V) 



less than Eq. 4.8. In this case, the slope of the fragment distribution matters. 
This scenario is represented as case B in Fig. 4.1. 

4.2.3 Fragment dominated regime 

The result obtained in the previous section may seem to be quite general. However, it 
does not hold for low .^-values as we will show in the following. 

In our case, the integrals do not diverge due to the finite integration bounds, however 
Makino et al. (1998) used and oo as lower and upper bounds for the integration and 
thus needed to investigate the convergence of the integral. They derived the following 
conditions: 

u-j-a + KO (4.18) 

7-a + 2>0. (4.19) 

where 7 gives the dependence of m2/mi in the /i-function of Eq. 4.20 (see Makino et al. 
1998): 



(4.20) 



The first condition (Eq. 4.18) considers the divergence towards the upper masses, 
whereas the second condition pertains the lower masses. We assume that Eq. 4.19 is 
satisfied and consider the case of decreasing ^ for a given by Eq. 4.17 which results in 
a steeper size distribution (where the mass is concentrated close to the upper end of 
the distribution). We see that for .^ < 1 + v — 27, Eq. 4.18 is no longer fulfilled. The 
behavior of the flux integral changes qualitatively: growth is no longer hierarchical, but 
it becomes dominated by contributions by the upper end of the integration bounds. 

Physically, this means that the total number of collisions of any grain is determined 
by the largest grains in the upper end of the distribution. Hence, smaller sized particles 
are predominantly refilled by fragmentation events of larger bodies instead of coagu- 
lation events from smaller bodies and they are predominantly removed by coagulation 
events with big particles (near the threshold rrif) instead of similar-sized particles. This 
corresponds to Case C in Fig. 4.1. 
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4.3 Simulation results interpreted 

To determine the resulting power-law distribution, we again focus on Eq. 4.5. The 
double integral of the flux F{m) can now be split into three separate integrals, 

F{m) = dnii \ dm2miCmi,m2n{mi)n{m2) 



dmi 

Jmo Jm—mi 



rm rmi 

+ dmi (lm2miCmi,m2'n{mi)n{m2) (4-21) 

J^ Jm—mi 

rm rnif 

+ dmi dm2miCmi^rn2'n{'mi)n{m2) 

according to whether m2 is larger or smaller than mi. 

It can be derived (see Appendix 4. A) that if the condition above (Eq. 4.18) is violated 
and if rrif » m, then the first and the third integral in Eq. 4.21 dominate the flux due 
to the integration until mf. In this cases, the flux F{m) is proportional to m?^'^^°'. 

A stationary state in the presence of fragmentation (Eq. 4.12) is reached if the fluxes 
cancel out, which leads to the condition 

a = C + 7. (4.22) 

This is the sweep-up regime where small particles are cleaned out by big ones (cf. Case 
C in Fig. 4.1). 



4.2.4 Summary of the regimes 

Summarizing these findings, we find that the resulting distribution is described by three 
scenarios (depicted in Fig. 4.1), depending on the slope of the fragment distribution: 

Case A (growth cascade): '^ > 2 a = ^^ 

Case B (intermediate regime): i^ — 27-|-l<^<2 a = '^^2 (4.23) 

Case C (fragment dominated) : ^<z^ — 27-I-I a = S^ + ^ 

4.3 Simulation results interpreted 

In this section, we will test the analytical predictions of the previous section by a co- 
agulation/fragmentation code (Birnstiel et al. 2010a, see also Brauer et al. 2008a). The 
code solves for the time evolution of the grain size distribution using an implicit integra- 
tion scheme. This enables us to find the steady-state distribution by using large time 
steps. In this way, the time evolution is not resolved, but the steady-state distribution 
is reliably and very quickly derived. 

We start out with the simplest case of a constant kernel and then - step by step - 
approach a more realistic scenario (in the context of a protoplanetary disk). In Sect. 4.4, 
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4. Dust size distributions in coagulation/fragmentation equilibrium 

we will then consider a kernel taking into account relative velocities of Brownian motion 
and turbulent velocities and also a fragmentation probability which depends on the 
masses and the relative velocities of the colliding particles. 

4.3.1 Constant kernel 

In the following section, we consider still the case of a constant kernel, but we will include 
fragmentation above particle sizes of 1 mm because this represents an instructive test 
case. 




10" 

grain size [cm] 



Figure 4.2: Grain size distributions^ for a constant kernel (i.e., z/ = 0) and different 
distributions of fragments. Tlie peak towards the upper end of the distribution is 
due to the fragmentation barrier (explanation in the text). The slope of the mass 
distribution corresponds to 6 — 3a. 



^It should be noted that in this work we will plot the distributions typically in terms of 
n{a) ■ m ■ a which is proportional to the distribution of mass. The advantage of plotting it 
this way instead of plotting n{a) is the following: when n{m) follows a power-law m^°' , then the 
grain size distribution n(o) describes a power-law with exponent 2 — 3a and the mass distribution 
exponent is 6 — 3a. For typical values of a, this is less steep and differences between a predicted 
and a real distribution are more prominent. 
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4.3 Simulation results interpreted 

We iteratively solve for the steady-state size distribution between coagulation and 
fragmentation. The outcome of these simulations for a constant kernel (i.e., 1^ = 0) are 
power-law distributions where the slope of the distribution depends on the fragmentation 
law (the slope ^, see Eq. 4.9). Figure 4.2 shows the corresponding size distributions for 
some of the different fragment distributions: the steepest distribution corresponds to the 
case of ^ = 0.5. For larger ^-values, the slope of the mass distribution flattens. In all 
cases, a bump develops towards the upper end of the distributions. The reason for this is 
the following: grains typically grow mostly through collisions with similar-sized or larger 
particles. Since the distribution is truncated at the upper end (defined as amax, see also 
Eq. 4.42), particles close to the upper end lack larger collision partners, the growth rate 
at these sizes would decrease if the distribution keeps its power-law nature. This, in 
turn means that the flux could not be constant below Omax- To keep a steady state, the 
number of particles at that point has to increase in order to replace the missing collision 
partners at larger sizes. 
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Figure 4.3: Exponent of grain size distributions for a constant kernel (i.e. z/ = 0) 
and different distributions of fragments (solid line) and the analytic solution for 
the growth cascade (Case A, dotted line), the intermediate regime (Case B, dashed 
line), and the fragment dominated regime (Case C, dash-dotted line). 



Figure 4.3 shows how the slope of the resulting distribution depends on the fragmen- 
tation slope ^, where the three previously discussed regimes can be identified: 
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4. Dust size distributions in coagulation/fragmentation equilibrium 

• Case A, growth cascade: as predicted, this scenario holds for values of ^ > 2 where 
most of the fragmenting mass is redistributed to fragments. This case corresponds 
to a "reversed" collisional cascade (just the direction in mass space is reversed as 
collisions mostly lead to growth instead of shattering). 

• Case B, the intermediate regime: when the fragment mass is more or less equally 
distributed over all sizes, mass gain by redistribution of fragments and the mass 
loss due to coagulation have to cancel each other. 

• Case C: most of the fragmented mass remains at large sizes. Therefore the mass 
distribution is dominated by the largest particles. In other words, growth is not 
hierarchical anymore. In this test case 7 = 0, and from Eq. 4.18 it follows that the 
transition between the intermediate and fragment-dominated regime lies at ^ = 1, 
which can be seen in Fig. 4.3. 

The measured slopes for m « nif are in excellent agreement with the model outlined in 
Sect. 4.2. 

4.3.2 Including settling effects 

The simplest addition to this model is settling: as grains become larger, they start to 
settle towards the mid-plane. However, turbulent mixing counteracts this systematic 
motion. The vertical distribution of dust in a settling-mixing equilibrium can (close to 
the mid-plane) be estimated by a Gaussian distribution with a size-dependent dust scale 
height H^. Smaller particles are well enough coupled to the gas to have the same scale 
height as the gas, Hp, while larger particles decouple and their scale height is decreasing 
with grain size, 

|^=./g, for St>at (4.24) 

lip Y ot 

(see e.g., Dubrulle et al. 1995; Schrapler &: Henning 2004; Youdin & Lithwick 2007) 
where the Stokes number St is a dimensionless quantity which describes the dynamic 
properties of a suspended particle. Very small particles have a small Stokes number and 
are therefore well coupled to the gas. Particles which have different properties (e.g., 
size or porosity) but the same St behave aerodynamically the same. In our prescription 
of turbulence, St is defined as the product of the particles stopping time Tst and the 
orbital frequency Q]^. We focus on the Epstein regime where the Stokes number can be 
approximated by 

St = f^k-rst-^5 (4.25) 

with Eg being the gas surface density and ps being the solid density of the particles 
which relates mass and size via m = Air/Spsa'^, 

Settling starts to play a role as soon as the Stokes number becomes larger than the 
turbulence parameter at which can be related to a certain size 

2a+ Tip- 
asett = -^^. 4.26 

vrps 
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4.3 Simulation results interpreted 

Eq. 4.26 only holds within about one gas pressure scale height because the dust vertical 
structure higher up in the disk deviates from the a Gaussian profile (see also DubruUe 
et al. 1995; Schrapler & Henning 2004; Duhemond & Dominik 2004). 

The mass-dependent dust scale height causes the number density distribution n{m) 
to depend on the vertical height, z. In disk-like configurations, it is therefore customary 
to consider the column density, 

rOO 

N{m) = n{m,z)dz. (4.27) 

J— 00 

Similar to Eq. 4.1, we can write the vertically integrated number of collisions as 

Crm,m2 ■ N{mi) ■ N{m2) dmi dm2, (4.28) 

which gives the total number of collisions that take place over the entire column of the 
disks. 

The dependence of the collisional probability on scale height, is now reflected in the 
modified kernel Cmi,m2 (see Section 3.1.2 for a derivation): 

^2vr {H! + HI) 

The point to realize here is that, due to the symmetry between Eq. 4.1 and Eq. 4.28, 
the analysis in Sect. 4.2 holds also for disk-like configuration, if the kernel is now replaced 
by Cmi,m2- The resulting exponent a then concerns the column density. 

If we consider the case that St > at and substitute Eq. 4.24 and Eq. 4.25 into 
Eq. 4.29, we find that 

^ _ ^m\,m2 

'^mi,m2 ~ 



2n {Hi + Hi) 



'^mi,m2 '"'1 " 



mi 



has an index v = 1/6 for grain sizes larger than Ogett and z/ = otherwise (it should 
be noted that Hi and H2 are the dust scale heights whereas h[m2/mi) represents the 
function defined in Eq. 4.2). 

The theory described in Sect. 4.2 is strictly speaking only valid for a constant u- 
index, but if this index is constant over a significant range of masses, then the local 
slope of the distribution will still adapt to this index. In the case of settling, we can 
therefore describe the distribution with two power-laws as can be seen in Fig. 4.4. 

The fact that the distribution will locally follow a power-law is an important require- 
ment for being able to construct fitting formulas which reproduce the simulated grain 
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4. Dust size distributions in coagulation/fragmentation equilibrium 

size distributions. It allows us in some cases to explain the simulation outcomes with 
the local kernel index (although coagulation and fragmentation are non-local processes 
in mass space, since each mass may interact with each other mass). A physically moti- 
vated recipe to fit the numerically derived distribution functions for the special case of 
^ = 11/6 is presented in Sect. 4.5. 
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Figure 4.4: Simulation result (solid line) and a "two power-law fit" to the vertically 
integrated dust distribution for a constant kernel with settling included. The 
dotted, vertical line denotes the grain size above which grains are affected by 
settling. 



4.3.3 Non-constant kernels 

We performed the same tests as in Sect. 4.3.1 also for non-constant kernels, i.e. kernels 
with 1/ > 0. The measured slopes for the cases of {u = 5/6, 7 = 0) (corresponding to the 
second turbulent regime, see Sect. 4.4.1), {u = 1/3,7 = 1/3), and {v = 1/6,7 = ~l/2) 
(i.e., Brownian motion, see Sect. 4.4.1) are shown in Fig. 4.5. Similar to Fig. 4.3, the 
distribution always follows the minimum index of the three different regimes (Cases 
A, B, and C). It should be noted that for all indices ^ of the fragmentation law, all 
processes - coagulation, fragmentation, and re-distribution of fragments - take place, 
but the relative importance of them is what determines the resulting slope. 
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4.4 Grain size distributions in circumstellar disks 

In Fig. 4.5, it can be seen that Case B (defined in Sect. 4.2.2) vanishes for the 
kernel in the upper panel, while it is present for a large range of S, for the kernel in the 
middle panel. This can be explained by the definitions of the three regimes which were 
summarized in Eq. 4.23: with (ly = 5/6,7 = 0) (cf- upper panel in Fig. 4.5), Case B is 
confined between ^ = 11/6 and 2. The grain size distribution, therefore, switches almost 
immediately from being growth dominated (Case A) to fragmentation dominated (Case 
C). In the case of a kernel with u = 1/3 and 7 = 1/3, this range spans from 2/3 to 2, as 
can be seen in the central panel of Fig. 4.5. 

The lower panel shows the distribution for a Brownian motion kernel (i.e., ly = 1/6 
and 7 = —1/2). The grey shaded area highlights the range of ^ values where our 
predictions do not apply, that is, Eq. 4.19 is not fulfilled any more. The resulting 
steady-state distributions are not power-law distributions anymore. 

4.4 Grain size distributions in circumstellar 
disks 

In this section, we will leave the previous "clean" kernels and focus on the grain size distri- 
bution in circumstellar disks including relative velocities due to Brownian and turbulent 
motion of the particles, settling effects, and a fragmentation probability as function of 
particle mass and impact velocity. 

Combining all these effects makes it impossible to find simple analytical solutions in 
the case of a coagulation-fragmentation equilibrium. We will therefore use a coagulation 
fragmentation code to find the steady-state solutions and to show how the steady-state 
distributions in circumstellar disks depend on our input parameters. 

We will first discuss the model "ingredients", i.e. the relative velocities and the 
prescription for fragmentation/sticking. Afterwards, we will define characteristic sizes 
at which the shape of the size distribution changes due to the underlying physics. In 
the last subsection, we will then show the simulation results and discuss the influence of 
different parameters. 

4.4.1 Relative velocities 

We will now include the effects of relative velocities due to Brownian motion, and due 
to turbulent mixing. The Brownian motion relative velocities are given by 

A.BM = >^-^(^^+^^^ (4.31) 

y vr nil '^'2 

where fcb is the Boltzmann constant and T the mid-plane temperature of the disk. 
Ormel &: Cuzzi (2007) have derived closed form expressions for particle collision velocities 
induced by turbulence. They also provided easy-to-use approximations for the different 
particle size regimes which we will use in the following. Small particles (i.e. stopping 
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Figure 4.5: Exponents of the grain size distributions for three different kernels as 
function of the fragment distribution index. The plot shows the simulation results 
(solid lines), the analytic solution for the growth cascade (Case A, dotted lines), the 
intermediate regime (Case B, dashed lines), and the fragment dominated regime 
(Case C, dash-dotted lines). The upper panel was calculated with a (z/ = 5/6,7 = 
0)-kernel (i.e. turbulent relative velocities, see Sect. 4.4.1), the middle panel with 
a (z/ = 1/3,7 = l/3)-kernel and the lower panel with a (z/ = 1/6,7 = — l/2)-kernel 
(i.e. Brownian motion relative velocities). In the grey shaded area, Eq. 4.19 is not 
fulfilled, and the simulation result is not a power-law like distribution. 
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Table 4.1: Kernel indices for the different regimes without settling. 



Regime 



Brownian motion regime 
Turbulent regime I 
Turbulent regime II 



7 upper end 



6 2 



Obt 
1 ai2 



a 



max 



time of particle « eddy crossing time) belong to the first regime of Ormel & Cuzzi (2007) 
with velocities proportional to 

Aiiioc|Sti-St2|. (4.32) 

The relative velocities of particles in the second turbulent regime of Ormel & Cuzzi 
(2007) are given by 

AniiocVStmax, (4.33) 

where Stmax is the larger of the particles Stokes numbers. Velocities in this regime show 
also a weak dependence on the ratio of the Stokes numbers which we will neglect in the 
following discussion. 

Together with the geometrical cross section cTgco = 7r(ai + 02)^, it is straight-forward 
to estimate the indices of the kernel, u and 7, as defined in Eq. 4.2 and 4.20 for all these 
regimes, without settling (assuming that only Brownian motion or turbulent motion 
dominates the relative velocities) . The indices for these three sources of relative particle 
motion are summarized in Table 4.1. If settling is to be included, then the u index for 
particle sizes above Ogctt has to be increased by i (see Sect. 4.3.2) while 7 remains the 
same. The i' and 7 indices for all three regimes can be found in Table 4.1. 

4.4.2 Fragmentation and cratering 

We introduce fragmentation and cratering according to the recipe of Birnstiel et al. 
(2010a): whether a collision leads to sticking or to fragmentation/cratering is determined 
by the fragmentation probability 

if An < U{ — 6u 

P{= ■{ I if Au > Uf (4.34) 

. 1-^^ else 

which means that all impacts with velocities above the critical break-up threshold Uf 
lead to fragmentation or cratering while impacts with velocities below Uf — 5u lead to 
sticking. The width of the linear transition region Su is chosen to be 0.2 Uf. 

The mass ratio of the two particles determines whether the impact completely frag- 
ments the larger body (masses within one order of magnitude) or if the smaller particle 
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Figure 4.6: Fiducial model and variations of the most important parameters: frag- 
mentation power-law index ,^, critical fragmentation velocity -Uf, turbulence param- 
eters at, surface density Sg, particle solid density ps, and mid-plane temperature 
T. The shape of the vertically integrated size distributions does not depend on 
the stellar mass or the distance to the star (only via the radial dependence of the 
parameters above). 
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4.4 Grain size distributions in circumstellar disks 

excavates mass from the larger body (masses differ by more than one order of magni- 
tude). 

In the case of fragmentation, the whole mass of both collision partners is redistributed 
to all masses smaller than the larger body according to Eq. 4.9. 

In the case of cratering, it is assumed that the smaller particle (with mass fnimp) 
excavates its own mass from the larger body. The mass of the impactor as well as the 
excavated mass is then redistributed to masses smaller than the impactor mass according 
to Eq. 4.9. Thus, the total redistributed mass equals 

n{m) ■ mdm = 2?7iimp (4.35) 

'mo 

and the mass of the larger body is reduced by 2 rriirap- 

Most parameters such as the fragmentation velocity or the amount of excavated 
material during cratering are not yet well enough constrained (for the most recent ex- 
perimental work, see Blum & Wurm 2008, Giittler et al. 2010, and references therein). 
Experiments suggest fragmentation velocities of a few m s^^ and fragment distributions 
with S, between 1 and 2 (Giittler et al. 2010 find ^ values between 1.07 and 1.37 for Si02 
grains). Simulations of silicate grain growth around 1 AU show that also bouncing (i.e. 
collisions without sticking or fragmentation) can play an important role (see Weidling 
et al. 2009; Zsom et al. 2010). However, changes in material composition such as organic 
or ice mantles or the monomer size are expected to change this picture. As there is still 
a large parameter space to be explored, we continue with the rather simple recipe of 
sticking, fragmentation and cratering outlined above. 

4.4.3 Regime boundaries 

From Eq. 4.23, we can calculate the slope of the distribution in the different regimes if 
we assume that the slope of the distribution at a given grain size always follows from 
the kernel index u. To construct a whole distribution consisting of several power- laws 
for each regime, we need to know where each of the different relative velocity regime 
applies. 

It is important to note that relative velocities due to Brownian motion decrease 
with particle size whereas the relative velocities induced by turbulent motion increase 
with particle size (up to St = 1). Therefore, Brownian motion dominates the relative 
velocities for small particles, while for larger particles, turbulence dominates. From 
numerical simulations, we found that at those sizes where the highest turbulent relative 
velocities (i.e. collisions with the smallest grains) start to exceed the smallest Brownian 
motion relative velocities (i.e. collisions with similar sized grains), the slope of the 
distribution starts to be determined by the turbulent kernel slope. By equating the 
approximate relative velocity of Ormel & Cuzzi (2007) and Eq. 4.31, the according grain 
size can be estimated to be 



obt 






1-2 



5 

(4.36) 
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4. Dust size distributions in coagulation/fragmentation equilibrium 

where we approximate the Reynolds number (i.e., the ratio of turbulent viscosity ft = 
aCgHp over molecular viscosity) near the disk mid-plane by 

Re^"^^^^. (4.37) 

2 finip 

Here, CHa is the cross section of molecular hydrogen (taken to be 2 x 10^^^ cm^) and 
H = 2.3 is the mean molecular weight in proton masses nip. 

Turbulent relative velocities strongly increase for grains with a stopping time that is 
larger or about the turn-over time of the smallest eddies. More specifically, the Stokes 
number of the particles at this change in the relative velocity is 

Sti2 = -^ = -Re-5, (4.38) 

Va tL Va 

where ti, = 1/$!^ and t^ = ^l ■ Re^2 are the turn-over times of the largest and the 
smallest eddies, respectively, and ilk is the Kepler frequency, ya is a factor which was 
approximated by Ormel & Cuzzi (2007) to be about 1.6. The corresponding grain size 
in the Epstein regime is therefore given by 

1 2S(T „ 1 , , 

ai2 = ^-Re^a, (4.39) 

Va vr/9s 

As mentioned above, the Brownian motion relative velocities of small grains decrease 
with their size. For larger sizes, the relative velocities due to turbulent motion are 
gaining importance, which are increasing with size until a Stokes number of unity. For 
typical values of the sound speed 

cs = .1^ (4.40) 

and the turbulence parameter at, the largest turbulent relative velocity Aumux ~ s/otCg 
exceed the critical collision velocity of the grains (which is in the order of a few m s^^) 
and therefore leads to fragmentation of the dust particles. In the case of very quiescent 
environments and/or larger critical collision velocities, particles do not experience this 
fragmentation barrier and can continue to grow. Hence, a steady state is never reached. 
The work presented here focuses on the former case where 

Alimax > Uf, (4.41) 

and grain growth is always limited by fragmentation. 

As relative turbulent velocities are (in our case) increasing with grain size, we can 
relate the maximum turbulent relative velocity and the critical collision velocity to derive 
the approximate maximum grain size which particles can reach (see Birnstiel et al. 2009) 

a^ax ^^^A. (4.42) 

iratPs ci 
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4.4 Grain size distributions in circumstellar disks 
4.4.4 Resulting steady-state distributions 

The parameter space is too large to even nearly discuss all possible outcomes of steady 
state grain size distributions. We will therefore focus on a few examples and rather 
explain the basic features and the most general results only. For this purpose, we will 
adopt a fiducial model and consider the influences of several parameters on the resulting 
grain size distribution: ^, Uf, at, Sg, Ps, and T (see Fig. 4.6). 

The flducial model (see the solid black line in Fig. 4.6) shows the following features: 
steep increase from the smaller sizes until a few tenth of a micrometer. This relates to 
the regime dominated by Brownian motion relative velocities. The upper end of this 
regime can be approximated by Eq. 4.36. The flatter part of the distribution is caused 
by a different kernel index v in the parts of the distribution which are dominated by 
turbulent relative velocities. The dip at about 60 /um (cf. Eq. 4.39) is due to the jump 
in relative velocities as the stopping time of particles above this size exceeds the shortest 
eddy turn-over time (see Ormel & Cuzzi 2007). 

The upper end of the distribution is approximately at Omax- The increased slope 
of the distribution and the bump close to the upper end are caused by two processes. 
Firstly, a boundary effect: grains mostly grow by collisions with similar or larger sized 
particles. Grains near the upper end of the distribution lack larger sized collision partners 
and therefore the number density needs to increase in order to keep the flux constant 
with mass (i.e. in order to keep a steady-state). Secondly, the bump is caused by 
cratering: impacts of small grains onto the largest grains do not cause growth or complete 
destruction of the larger bodies, instead they erode them. Growth of these larger bodies 
is therefore slowed down and, similar to the former case, the mass distribution needs to 
increase in order to fulfill the steady-state criterion ("pile- up effect"). 

The upper left panel in Fig. 4.6 shows the influence of the distribution of fragments 
after a collision event: larger values of ^ mean that more of the fragmented mass is redis- 
tributed to smaller sizes. Consequently, the mass distribution at smaller sizes increases 
relative to the values of smaller ^ values. 

The strong influence of the fragmentation threshold velocity Uf can be seen in the 
upper right panel in Fig. 4.6: according to Eq. 4.42, an order of magnitude higher Uf 
leads to a 100 times larger maximum grain size. 

The grain size distributions for different levels of turbulence are shown in the middle 
left panel of Fig. 4.6. The effects are two-fold: flrstly, an increased at leads to increased 
turbulent relative velocities, thus, moving the fragmentation barrier amax to smaller sizes 
(cf. Eq. 4.42). Secondly, a larger at shifts the transition within the turbulent regime, 
ai2, to smaller sizes. Consequently, the second turbulent regime gains importance as at 
is increased since its upper end lower boundary extend ever further. 

The middle right panel in Fig. 4.6 displays the influence of an decreased gas surface 
density Sg (assuming a fixed dust-to-gas ratio). It can be seen that not only the to- 
tal mass is decreased due to the fixed dust-to-gas ratio but also the upper size of the 
distribution amax decreases. This is due to the coupling of the dust to the gas: with 
larger gas surface density, the dust is better coupled to the gas. This is described by a 
decreased Stokes number (see Eq. 4.25) which, in turn, leads to smaller relative velocities 
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4. Dust size distributions in coagulation/fragmentation equilibrium 



and hence a larger amax- Interestingly, the shape of the grain size distribution does not 
depend on the total dust mass, but on the total gas mass, as long as gas is dynamically 
dominating (i.e., Sg » S^). If more dust were to be present, grains would collide more 
often, thus, a steady-state would be reached faster and the new size distribution would 
be a scaled-up version of the former one. However the velocities at which grains collide 
are determined by the properties of the underlying gas disk. In this way, the dust grain 
size distribution is not only a measure of the dust properties, but also a measure of the 
gas disk physics like the gas density and the amount of turbulence. 

The shape of the size distribution for different grain volume densities ps does not 
change significantly. However most regime boundaries (osett, ai2, and amax) are inversely 
proportional to ps (because of the coupling to the gas, described by the Stokes number). 
A decrease (increase) in p^ therefore shifts the whole distribution to larger (smaller) sizes 
as can be seen in the lower left panel in Fig. 4.6. 

The upper end of the distribution, amax, is inversely proportional to the mid-plane 
temperature T (as in the case of the turbulence parameter) whereas the transition be- 
tween the different turbulent regimes ai2 does not. Therefore, increasing the temperature 
decreases Omax in the same way as decreasing Sg does. However ai2 is not influenced by 
temperature changes, therefore the shape of the size distribution changes in a different 
way than in the case of changing Eg as can be seen by comparing the middle right and 
the lower right panels of Fig. 4.6. 



4.5 Fitting formula for steady-state 
distributions 

In this section, we will describe a simple recipe which allows us to construct vertically 
integrated grain size distributions which fit reasonably well to the simulation results 
presented in the previous section. 

The recipe does not directly depend on the radial distance to the star or on the 
stellar mass. A radial dependence only enters via radial changes of the input parameters 
listed in Table 4.3. This recipe has been tested for a large grid of parameter values^, 
however there are some restrictions. 

4.5.1 Limitations 

These fits strictly apply only for the case of ^ = 11/6. In this case, the slopes of the 
distribution agree well with the predictions of the intermediate regime (Case B, defined 
in Sec. 4.2.2). For smaller values of ^, the slopes do not strictly follow the analytical 
predictions. This is due to the fact that we include cratering, which is not covered by 
our theory. Erosion is therefore an important mode of fragmentation: it dominates over 
complete disruption through the high number of small particles (see also Kobayashi &: 



^See www. mpia.de/distribution- fits 
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10^ 
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10^ 
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fit: Ait^i 
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Figure 4.7: Comparison of the turbulent relative velocities of Ormel & Cuzzi 
(2007) to the fitting formula used in the recipe (see Eqs. 4.47 and 4.44) for grain 
size distributions. The largest error in the resulting upper grain size ap derived 
from the fitting formula is about 25%. 



Tanaka 2010) and it is able to redistribute significant amounts of mass to the smallest 
particle sizes. 

One important restriction for this recipe is the upper size of the particles Omax: it 
needs to obey the condition of Eq. 4.41, since otherwise, particles will not experience 
the fragmenting high-velocity impacts and a steady state will never be reached since 
particles can grow unhindered over the meter-size barrier. 

There are also restrictions to a very small amax^ if ^max is close to or even smaller 
than ai2, then the fit will not represent the true simulation outcome very well. In this 
case, the upper end of the distribution, and in more extreme cases the whole distribution 
will look much different. Thus, the sizes should obey the condition 



5 fim < ai2 < «n 



(4.43) 



where each inequality should be within a factor of a few. 
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Figure 4.8: Step-by-step construction of the fit distribution for the following pa- 
rameters: Eg = 20 g cm^^, Sd = 0.2 g cm^^, at = 1 x. 10^'', -Uf = 1 m s^^, 
^ = 1.833, T = 50 K, ps = 1-6 g cm^^. The upper left panel shows the distribution 
after step 5: each interval obeys a different power-law, the distribution is continu- 
ous apart from a jump at ai2- The upper right panel displays the fit including and 
increase at the upper end, according to step 6. The bump caused by cratering (cf. 
Eq. 4.52) is shown in the bottom left panel while the bottom right panel compares 
the final fit distribution (solid curve) to the simulation result (dashed curve). The 
vertical lines correspond to the regime boundaries: Obt (solid), Cgctt (dash-dotted), 
ai2 (dashed) and ap (dotted). 
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4.5.2 Recipe 



The following recipe calculates the vertically integrated mass distribution of dust grains 
in a turbulent circumstellar disk within the the above mentioned limitations. The recipe 
should be applied on a logarithmic grain size grid a^ with a lower size limit of 0.025 nm 
and a fine enough size resolution (oi+i/aj < 1.12). For convenience, all variables are 
summarized in Table 4.3. The steps to be performed are as follows: 

1. Calculate the grain sizes which represent the regime boundaries aeT, 012, Osctt 
which are given by Eqs. 4.36, 4.39, and 4.26. 

2. Calculate the turbulent relative velocities for each grain size. For this, we approx- 
imate the equations which are given by Ormel & Cuzzi (2007). Collision velocities 
with monomers are given by 



Auf 



^gas 



Re4 • (Sti - Sto) 

(1 - e) • Re3 • (Sti 
+ e ■^/3^Sti 



for Oj < ai2 



Sto) for ai2 < Cj < 5 ai2 



for Oj > 5 ai2 



(4.44) 



where Re is the Reynolds number (see Eq. 4.37), and Stj and Sto are the Stokes 
numbers of Oj and monomers, respectively (cf. Eq. 4.25). Ugas is given by 



u, 



gas 



:"t, 



and the interpolation parameter e is defined as 

_ Oj — ai2 
4ai2 

and collisions with similar sized bodies are approximated as 

for Oj < ai2 



(4.45) 



(4.46) 



Anr 



(4.47) 



i • Au™°" 



for Oj > 012 



A comparison between these approximations and the formulas of Ormel &: Cuzzi 
(2007) is shown in Fig. 4.7. 

3. Using Eqs. 4.44 and 4.47 for the relative velocities, and the transition width 5u = 
0.2uf, find the grain sizes which correspond to the following conditions: 

• cl: particles above this size experience impacts with monomers with veloci- 
ties of Auf^°^ ^ U{ — 6u (i.e., cratering starts to become important). 



119 



4. Dust size distributions in coagulation/fragmentation equilibrium 

• ap: particles above this size experience impacts with equal sized grains with 
velocities of Au'^'^ ^ Ui — du (i.e., fragmentation becomes important). 

• Or: particles above this size experience impacts with similar sized grains with 
velocities of An^^ ^ Uf (i.e. every impact causes fragmentation/cratering) . 

4. Calculate the factor J according to the recipe 

V = cJ'^!^^^)'.jy^ (4.48) 

V at CTHa / Y 4 Eg ya 

1 
J= I 2.5-^+ I 1.1«+ (1 + 2V3-) 1 I , (4.49) 



where fj, = 2.3, a^^ = 2 x 10 ^^ cm^ and ya = 1.6. 



5. The power-law indices of the mass distribution 6i for each interval between the 
regime boundaries (asT, ai2, Osett, op) are calculated according to the intermediate 
regime (cf. Eq. 4.17). The slopes have to be chosen according to the kernel regime 
(Brownian motion, turbulent regime 1 or 2) and according to whether the regime is 
influenced by settling or not (i.e., if a is larger or smaller than agett)- The resulting 
slopes of the mass distribution are given in Table 4.2. The first version of the fit 
/(oj) (where Oj denotes the numerical grid point of the particles size array) is now 
constructed by using power-laws (<xa^') in between each of the regimes, up to ap. 
The fit should be continuous at all regime boundaries except a drop of 1/J at the 
transition at ai2- An example of this first version is shown in the top left panel of 
Fig. 4.8. 

6. Mimic the cut-off effects which cause an increase in the distribution function close 
to the upper end: linearly increase the distribution function for all sizes ainc < 
a < ap: 

f{ai) - /(ai) • (2 - °' ~ "M , (4.50) 

\ Oinc — ap / 

with 

Oinc = 0.3 ap. (4.51) 

The resulting fit after this step is shown in the upper right panel of Fig. 4.8. 

7. The bump due to cratering is mimicked by a Gaussian, 

6(a.) = 2 • /(aL) • exp ( -i^i^^ J , (4.52) 



where a is defined as 



min(|aR-ap|,|aL-ap|) 
a = , , (4.5d) 
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Table 4.2: Power-law exponents of the distribution n{m) ■m-a. The slopes were cal- 
culated using the formula for a coagulation/fragmentation equilibrium (Eq. 4.17). 
Within each regime of relative velocities, it has to be differentiated whether grains 
are influenced by settling or not. 



Regime 


6 


i 


di < Osctt 


O-i > flsctt 


Brownian motion regime 


3 

2 


5 
4 


Turbulent regime I 


1 
4 





Turbulent regime II 


1 
2 


1 

4 



but should be limited to be at least 



a > 0.1 Op. 



(4.54) 



The fit T{ai) is now constructed by using the maximum of /(aj) and 6(aj) in the 
following way: 



J"(aj) = < 



f{ai) if Oj ig QL 

max(/(aj),6(ai)) if cl ^ Oj ^ ap 

b{ai) if ap ^ ai ^ or 

else 



(4.55) 



9. Finally, the fit needs to be normalized to the dust surface density at the given 
radius. J-" is a (yet un- normalized) vertically integrated mass distribution (shown 
in the bottom right panel of Fig. 4.8). To translate this mass distribution to a 
vertically integrated number density distribution N{a), we need to normalize it as 



N{a) 



S, 



Ln -^(") dl^« 



J-(a) 



m ■ a 



(4.56) 



4.6 Conclusions 

In this work, we generalize the analytical findings of previous works to the case of 
grain size distributions in a coagulation/fragmentation equilibrium. Under the assump- 
tion that all grains above a certain size, Omax) fragment into a power-law distribution 
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Table 4.3: Definition of tlie variables used in this chapter. The variables grouped 
as "input variables" are the parameters of the fitting recipe in Sect. 4.5. "Other 
variables" summarizes the definitions of all other variables used in this work. 





variable 


definition 


unit 


s 


at 


Turbulence strength parameter 


- 


1 


U( 


Fragmentation threshold velocity 


cm s"""^ 


u 

g 


Sg 


Gas surface density 


g cm"^ 




Ed 


Dust surface density 


g cm"^ 


? 


Power-law index of the mass distribution of fragments, see Eq. 4.9 


- 


1 — 1 


T 


Mid-plane temperature 


K 




Ps 


Volume density of a dust particle 


g cm-^* 




a 


Slope of the number density distribution n(7n)ocm~™, see Eq. 4.3 


- 




IBT 


Eq. 4.36 


cm 




ai2 


Eq. 4.39 


cm 




^max 


Eq. 4.42 


cm 




aL 


Left boundary of the bump function fe(a), see Sect. 4.5, paragraph 3 


cm 




ap 


Peak size of the bump function &(a), see Sect. 4.5, paragraph 3 


cm 




an 


Right boundary of the bump function &(a), see Sect. 4.5, paragraph 3 


cm 




<^sctt 


Eq. 4.26 


cm 




^inc 


Eq. 4.51 


cm 




b{a) 


Bump function, see Eq. 4.52 


- 




^mi ,1712 


Collision kernel, see Eq. 4.2 


cm^ s"^ 




^rrii ,m2 


Collision kernel including settling effects, see Eq. 4.2 


2 —1 

cm s 


m 

^ 


Cs 


sound speed, see Eq. 4.40 


cm s~^ 


Xi 
.s 


Si 


Slopes of the fit-function, see Table 4.2 


- 


'S 


5u 


Width of the transition between sticking and fragmentation, taken to 


cm s""*^ 


tH 




be 0.2 Ui 






7 


Power-law index of the function h{m2/mi) for large ratios of m2/mi, as 


- 


o 




defined in Eq. 4.20 






K 


Integral defined in Eq. 4.6 


- 




feb 


Boltzmann constant 


erg K~^ 




mo 


Monomer mass 


g 




nil 


Mass above which particles fragment 


g 




M 


Mean molecular weight in proton masses, taken to be 2.3 


- 




rrip 


Proton mass 


g 




n{rn) 


Number density distribution 


g cm"^ 




N{m) 


Vertically integrated number density distribution 


g cm"^ 




V 


Degree of homogeneity of the kernel as defined in Eq. 4.2 


- 




Re 


Reynolds number, see Eq. 4.37 


- 




St 


Mid-plane Stokes number in the Epstein regime, see Eq. 4.25 


- 




0"H2 


Cross-section of molecular hydrogen 


2 

cm 




Aur°" 


Relative velocities between monomers and grains of size at, see Eq. 4.44 


cm s"""^ 




A^r 


Relative velocities between grains of size Ui, see Eq. 4.47 


cm s"""^ 
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4. A Derivation of the fragment dominated size distribution 

nf(r7T,)ocr7T,^^, we derived analytical steady-state solutions for self-similar kernels and de- 
termined three different cases (see 4.2.4). 

Results show that dust size distributions in circumstellar disks do not necessarily 
follow the often adopted MRN power-law distribution of n(a)oca^^'^ (see Mathis et al. 
1977; Dohnanyi 1969; Tanaka et al. 1996; Makino et al. 1998; Garaud 2007) when both 
coagulation and fragmentation events operate. We performed detailed simulations of 
grain growth and fragmentation to test the analytical predictions and found very good 
agreement between the theory and the simulation results. 

We applied the theory to the gaseous environments of circumstellar disks. Unlike 
the models of Garaud (2007), the upper end of the size distribution is not limited by 
the growth timescale but by fragmentation as relative velocities increase with grain size 
and reach values large enough to fragment grains. The shape of the dust distribution 
is determined by the gaseous environment (e.g., gas surface density, level of turbulence, 
temperature and others) since the gas is dynamically dominating as long as the gas 
surface density significantly exceeds the dust surface density. The total dust mass merely 
provides the normalization of the distribution and the timescale in which an equilibrium 
is reached. The results presented in this work show that the physics of growth and 
fragmentation directly link the upper and the lower end of the dust distribution in 
circumstellar disks. 

A ready-to-use recipe for deriving vertically integrated dust size distributions in 
circumstellar disks for a fixed value of ^ = 11/6 is presented in Sect. 4.5. Although the 
collision kernel in circumstellar disks is complicated, we found good agreement with our 
fitting recipe for a fragment distributions with ^ = 11/6. The recipe can readily be used 
for further modeling such as disk chemistry or radiative transfer calculations. 

4. A Derivation of the fragment dominated size 
distribution 

In this section, we will derive the slope of the size distribution which is dominated by the 
largest particles. Since in our scenario, the integrals are confined between the monomer 
mass rriQ and the largest particles at the fragmentation barrier rrif, the integrals do not 
diverge as in the scenario of Tanaka et al. (1996) and Makino et al. (1998), who consider 
integration bounds of zero and infinity. However, if the first condition of Makino et al. 
(1998) is not fulfilled, then the mass fiux (cf. Eq. 4.5) is dominated by the upper bound 
of the integral. This is the case which we will consider in the following. 

As noted in Sect. 4.2.3, the fiux integral (Eq. 4.5) can be split into three separate 
integrals, in order to distinguish cases of m2 > rrii or m2 < mi, 

F{m) ^h + l2 + h, (4.57) 

where /i, I2, and I3 correspond (from left to right) to the three integrals defined in 
Eq. 4.21. We will now evaluate these integrals in the limits of rriQ « m « rrif, using the 
limiting behavior of h{m2/rai) as given in Eq. 4.20. 
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4.A.1 First integral 

Carrying out the integration over 777-2 in -^1 leads to 

A^ ■ ho r™/2 



/l 



u — J — a + 



dnii mj 

mo 



7—0+1 



(4.58) 



\m 



V — c— a+1 



(m — 777l) 



J/— 7— a+1 



from which we can derive Eq. 4.18, the first convergence criterion of Makino et al. (1998). 
We consider the case where this condition does not hold, i.e.. 



jv — 7 — a + l>0. 



(4.59) 



Then, the TTif term in brackets dominates over the other term. Thus, the term in brackets 
is constant and carrying out the integration yields 



A^-ho- mp^-"+^ 



/rr7\7-a+2 2' 



(z^ - 7 - a + 1) • (7 - Q + 2) 
Now, if the second condition, Eq. 4.19 holds, using mg « m, we derive 

^2 . /jn . ,„J'-T-°+i /rr7\7-"+2 



h 



iQ • 777^ 



(z/ - 7 - a + 1) • (7 - a + 2) 



(4.60) 



(4.61) 



4. A. 2 Second integral 

We rewrite I2 using the dimensionless variables x\ = 7771/777 and X2 = 7712/771 which yields 

l2 = A'^-ho- 7773+'^-2" r dxi r dx2 x}+'' 

By integrating over X2, we derive 



hu-'j-a 7-Q 

* Xo 



(4.62) 



A'^-hn- m3+^-2° n 



7 — a + 



- — I dxix}+"-^-" • [xr"+^ - (1 



Xl 



>7— a+1 



(4.63) 



:=D 



The term in square brackets can be split into a sum of integrals, the first of which is 
straight-forward to evaluate as 



D 



l\3+u-2a 



u-2a + 3 



c..^ 



+ u — ^ — a 



•(1-xi) 



7 — Q + 1 



(4.64) 



while the second can be identified as a sum of a Beta function B[a, b) and an incomplete 
Beta function Bi (a, 6), 



■l\a+fe-l 






(4.65) 
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The condhions from above, 

a = 2 + i/-7-a >0 (4.66) 

6 = 2 + 7-a >0 (4.67) 

assure that the Beta functions, and thus I2, are reaL The numerical value of I2/D for a 
range of values of a and h are shown in Fig. 4.9. 

I2/B 




5 b 

Figure 4.9: Integral I2/D as function of a and h. 



4. A. 3 Third integral 

Similarly, we can rewrite /a as 



h = A'- hn 



nm nrrii 

rm 

— ■ dmi m^ 
1 Jr 



J 1+7—0 J/— 7—0 



■ rrio 



A' -ho 

V — ^ — a + 

' 2 

If Eq. 4.66 holds, then from mf » m follows 
^ ~ (z^ - 7 - Q + 1) • (7 - a + 2) ' 



1+7-a _ { ^^-7-«+i _ /m\^-7-"+i 



7-0+2^ 



• m 



7-a+2 



(4.68) 



(4.69) 
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4. A. 4 Deriving the steady-state distribution 

The first and the third integrand /i and I^ show the same mass dependence and can be 
summed up to 

^2^ocmp^-"+^ • m^-"+2 (4.70) 
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while I2 is proportional to 



^^ ocm^+'^-^a 



A^-h 







i^— a— 7+1 



oc I — 1 rrif m' ^ 



(4.71) 



and therefore 



/2 fmY-''-''+' , , 

-^ — • (4-72) 



I1 + /3 V"^f/ 

Since the constants of proportionality are factors of order unity and ra < rrif, the integrals 
/i + /3 are much larger than I2, therefore, the flux F{m) is proportional to rrCi^"^^"^. 

In the case of a steady state, this flux and the downward flux of fragments (which is 
proportional to m?^^) have to cancel each other. Therefore, the exponents of the mass 
dependence need to cancel out, i.e. 

« = 7 + e, (4.73) 

which is the slope of the steady-state condition in the fragmentation dominated regime 
(Case C in Fig. 4.1). 



CHAPTER 

Dust retention in protoplanetary 

disks 

Published as Birnstiel, Dullemond, & Brauer 2009, A&A, 503, L5 

Abstract 

Protoplanetary disks are observed to remain dust-rich for up to several 
million years. Theoretical modeling, on the other hand, raises several ques- 
tions. Firstly, dust coagulation occurs so rapidly, that if the small dust 
grains are not replenished by collisional fragmentation of dust aggregates, 
most disks should be observed to be dust poor, which is not the case. Sec- 
ondly, if dust aggregates grow to sizes of the order of centimeters to meters, 
they drift so fast inwards, that they are quickly lost. We attempt to verify if 
collisional fragmentation of dust aggregates is effective enough to keep disks 
'dusty' by replenishing the population of small grains and by preventing 
excessive radial drift. With a new and sophisticated implicitly integrated 
coagulation and fragmentation modeling code, we solve the combined prob- 
lem of coagulation, fragmentation, turbulent mixing and radial drift and at 
the same time solve for the 1-D viscous gas disk evolution. We find that for 
a critical collision velocity of 1 m/s, as suggested by laboratory experiments, 
the fragmentation is so effective, that at all times the dust is in the form 
of relatively small particles. This means that radial drift is small and that 
large amounts of small dust particles remain present for a few million years, 
as observed. For a critical velocity of 10 m/s, we find that particles grow 
about two orders of magnitude larger, which leads again to significant dust 
loss since larger particles are more strongly affected by radial drift. 
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Dust retention in protoplanetary disks 
5.1 Introduction 

The quest for a comprehensive understanding of the formation of planets in the dusty 
circumsteUar disks around many young stars has gained significant impetus. This is 
partly because of developments in the numerical modeling of these processes, but also 
because of the high-quality infrared and (sub-)millimeter data that are now available of 
hundreds of T Tauri and Herbig Ae star disks from observatories such as the Spitzer 
Space Telescope (e.g. Furlan et al. 2006; Kessler-Silacci et al. 2006), the Very Large Tele- 
scope (e.g. van Boekel et al. 2004), the Submillimeter Array (e.g. Andrews & Williams 
2007) and the Plateau de Bure Interferometer (e.g. Pietu et al. 2007). The challenge for 
theoretical astrophysicists is now not only to come up with a plausible model of how to 
overcome the various problems in our current understanding of how planets form, but 
also to ensure that new models are consistent with the observational data obtained from 
these planetary birthplaces. 

One of the challenges of theories of planet formation is what is often called the 
'meter-size barrier'. It was shown by Weidenschilling (1977) that as dust particles grow 
by coagulation, they acquire increasingly large relative velocities with respect to other 
particles in their vicinity, as well as a systematic inward drift velocity as originally 
described by Whipple (1972). The high relative velocities lead to destructive collisions 
(Blum &; Wurm 2008) , and thereby limit the growth of the aggregates to some maximum 
size. The radial drift also leads to a loss of solids toward the star well before large bodies 
can even form (Brauer et al. 2008a, hereafter B08a). 

Dullemond & Dominik (2005) showed that if, hypothetically, coagulation can proceed 
without limit, i.e., if no fragmentation occurs and no radial drift is present, the growth is 
so efficient that within 10^ years very few small (a < 1mm) opacity bearing grains remain 
in the disk. The disk becomes optically thin and the mid- to far-infrared fiux of the disk 
drops to levels well below that observed in the majority of gas-rich circumsteUar disks 
around pre-main sequence stars. By excluding other mechanisms, Dominik &: Dullemond 
(2008) argue that fragmentation is the most promising process for maintaining a high 
abundance of fine-grained dust in evolved disks. 

However, it has never been demonstrated explicitly that fragmentation is efficient 
enough to retain the population of small grains at the required levels. Moreover, if 
we include the effect of radial drift, it remains to be seen whether this radial drift is 
suppressed sufficiently by keeping the grains small through fragmentation. It is the goal 
of this chapter to investigate this issue and to determine how efficient fragmentation 
should be to agree with observations. 



5.2 Model 

The model presented here is a combination of a ID viscous gas disk evolution code 
and a dust evolution code treating the turbulent mixing, radial drift, coagulation and 
fragmentation of the dust. A detailed description of the model can be found in Chapters 2 
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5.2 Model 



and 3. Although the evolution of dust can influence the gas disk evolution, we neglect 
this eff'ect in our model. 

The gas disk model is similar to that described in Hueso &: Guillot (2005). The 
mid-plane temperature of the disk is approximated by following Nakamoto &: Nakagawa 
(1994), taking irradiation by the central star and viscous heating into account. The 
viscous evolution of the disk surface density Sg is described by (see Lynden-Bell & 
Pringle 1974), 

^_l^(v rn.)=0, (5.1) 

where the gas radial velocity Ug is given by 

ng = -^-^-(Sg.gVr^). (5.2) 

fg is the gas turbulent viscosity, 

Ug = aCsHp, (5.3) 

Cs denotes the sound speed, Hp the pressure scale height, and a is the turbulence pa- 
rameter (Shakura &: Sunyaev 1973). 

Grains are affected by a systematic radial inward drift related to headwind caused 
by the pressure-supported gas (Weidenschilling 1977), 

where St is the Stokes number of the particle (a dimensionless representation of the 
particle size), Pg is the gas pressure, pg is the gas volume density, and O^ is the Kepler 
frequency. In the Epstein regime, St is given by 

St = ^", (5.5) 

Sg2' ^ ^ 

where ps is the solid density of the dust grains. 

The second contribution to radial dust velocities is the gas drag due to the radial 
motion of the gas, 

Wdrag = g^, (5.6) 

where Sc = 1 -I- St is the Schmidt number, following Youdin & Lithwick (2007). 

The coupling to the gas turbulence leads to turbulent mixing of each dust species 
with a diffusion constant that is taken to be 

The dust grain number density n{m,r,z) is a function of mass m, radius r, height 
above the mid-plane z, and time. If we define 

poo 

a{'m,r) = n{'m,r, z)m'^dz, (5-8) 

J— 00 
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Dust retention in protoplanetary disks 

the vertically integrated time-evolution of this distribution can now be described by a 
general two-body process and an advection-diffusion equation, 

'^ g^'^ = JJq K{m,m',m") a{m',r) a{m",r)dm' dm" 

-liri^' '^im,r) ■ (Udrag + ^tdrift)] (5.9) 

Here, the right-hand side terms (from top to bottom) correspond to coagulation/frag- 
mentation, advection, and turbulent mixing. 

Since the detailed version of the coagulation/fragmentation equation is lengthy (Wei- 
denschilling 1980; Nakagawa et al. 1981; Dullemond & Dominik 2005; B08a), we include 
here only the general integral of a two-body process (first term on the right hand side). 
For a detailed description of the physics of coagulation/fragmentation, we refer to B08a; 
the numerical implementation is discussed in Chapter 3. 

The physical effects that produce the relative particle velocities considered here are 
Brownian motion, relative radial motion, vertical settling (see B08a), and turbulent 
motion (see Ormel & Cuzzi 2007). 

If particles collide with a relative velocity higher than the critical velocity Uf, they 
either fragment into a power-law size distribution of fragments (i.e., n{m)ccm^^'^^, see 
B08a) or the smaller body excavates mass from the larger one (cratering). We assume 
that the amount of excavated mass equals the mass of the smaller body. The fragmenta- 
tion velocity U{ is a free parameter in our model, and unless otherwise noted is taken to 
be 1 m/s, as suggested by experiments (e.g., Blum & Muench 1993) and by theoretical 
modeling (e.g., Stewart & Leinhardt 2009). 

We note that the velocities and therefore also the relative velocities are dependent 
on the particle size. If we assume that relative particle velocities are dominated by 
turbulent motion, 

AuturbOCVoSt Cg, (5.10) 

then we can estimate the Stokes number of the largest particles by equating the relative 
velocity and fragmentation velocity. 

With Eq. 5.5, we can relate this maximum 'dimensionless size' Stmax to a particle size, 

^ - ^ (5.12) 



"■max — 9 • 

Traps Cg 

It should be noted that Eq. 5.10 is an approximation of the equations derived in Ormel &: 
Cuzzi (2007). The detailed relative velocities of particles depend on the Stokes numbers 
of both particles. Hence, Eq. 5.12 is an order of magnitude estimate. 
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5.3 Results 



The combined equations for coagulation, fragmentation, and radial motion due to 
drift, drag, and mixing are solved using the technique of implicit integration (similar to 
B08a). We use the SixPACK F-90 linear algebra package^ to solve the sparse matrix 
equation. 



5.3 Results 

The simulation results presented here begin with a 0.07 Mq disk of power-law surface 
density S(r)ocr^^ (from 0.05 to 150 AU) around a 0.5 Mq star. Unless otherwise noted, 
we used a turbulence parameter of 10^^. The three models differ in the prescription 
of fragmentation. In model 1, we accounted only for growth without fragmentation. In 
model 2, particles fragment at Uf = 10 m/s, whereas in model 3, the fragmentation speed 
is taken to be 1 m/s. Selected results of all three models are depicted in Figure 5.1. 

The top row shows the evolution in the surface density distribution of dust because 
of coagulation, radial mixing, and radial drift in an evolving protoplanetary disk, cor- 
responding to model 1. The dashed line denotes the grain evaporation radius, which 
moves inwards as the temperature falls (lower surface density caused by accretion pro- 
duces less viscous heating) . The evaporation radius in our model is defined as the radius 
where the temperature rises above 1500 K assuming that all particles evaporate at this 
temperature. The solid line denotes the dust grain size, which corresponds to a Stokes 
number of unity. It can be seen that particles grow to decimeter size and then quickly 
drift inside the evaporation radius. 

The results of simulations that also include grain fragmentation are plotted in the 
second and third row of Figure 5.1. In the innermost regions of model 2, relative particle 
velocities are high enough (as a result of shorter dynamical times) to fragment particles 
already at relatively small Stokes number. However, the maximum size of the grain 
distribution is still strongly affected by radial drift, especially at larger radii. Larger 
particles are subject to radial drift and therefore drift to smaller radii, where they are 
pulverized. 

In contrast, fragmentation becomes the main limiting factor for particle growth 
throughout the disk if the critical fragmentation velocity is 1 m/s (model 3). The result 
is that particles fragment at a small Stokes number. Particles with a Stokes number less 
than about 10^^ are still strongly coupled to the gas and not subject to strong radial 
drift (see Brauer et al. 2007). Therefore, significant amounts of dust can be retained for 
several million years if fragmentation is included in the calculations. 

Figure 5.2 compares the evolution in the dust-to-gas ratio (solid line) and the 0.55 
yum optical depth at 10 AU (dashed line) of all three simulations. The optical depth is 
defined as 



/■oo 

a{m) Ky{m)dm. (5.13) 

Jo 



^available from www . engineers . auckland . ac . nz/~snor007 
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Figure 5.1: Snapshots of the vertically integrated dust density distribution of 
model 1 (only coagulation), model 2 (fragmentation at 10 m/s) and model 3 (frag- 
mentation at 1 m/s). The grain size is given by a = (3?7i/4'7rps)^/^. The dashed 
line shows the evaporation radius within which no coagulation is calculated. The 
solid line shows the particle size corresponding to a Stokes number of unity. 
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5.4 Discussion and conclusions 

It can be seen that the dust-to-gas ratio dechnes draniaticaUy in the first two cases. The 
optical depth in model 1 drops on even shorter timescales than the dust-to-gas ratio 
because not only the dust loss but also the dust growth causes the optical depth to 
decrease. 

Model 2 can retain the optical depth longer because small particles remain (due to 
fragmentation), providing enough surface to keep the disk optically thick for about 0.9 
Myr. The initial dip in the optical depth is caused by the initial condition: since initially, 
all mass is in small grains, the optical depth decreases as particles grow to larger sizes 
until fragmentation sets in, which increases the optical depth until a semi-equilibrium 
of growth and fragmentation is reached. 

If particles are already fragmented at 1 m/s (model 3), the dust-to-gas ratio declines 
only slightly, after 2 My, only about 4% of the dust being lost (compared to 96% and 94% 
in models 1 and 2). The optical depth decreases on longer timescales than in model 2 
since the dust mass is lost only on the viscous timescale. 

This mechanism of grain retention relies on the particles remaining small due to 
fragmentation. From Eq. 5.11 and 5.12, it follows that the maximum size or Stokes 
number depends strongly on the critical velocity Uf and less strongly on a. If ui is 10 
m/s instead of 1 m/s, the corresponding particle size is 100 times larger meaning that the 
largest particles are now also experiencing the radial drift barrier and therefore drifting 
quickly towards the central star. 

5.4 Discussion and conclusions 

We have shown that grain fragmentation plays an important role, not only in keeping 
the disk optically thick (by replenishing small grains) but also in maintaining a relatively 
high abundance of solids of all sizes. If particles fragment at relatively small sizes, they 
can stay below the size-regime where they are affected by radial drift. 

The simulation results presented here suggest that fragmentation of relatively low 
critical velocity (a few m/s) is needed to retain the dust required for planet formation, 
although the formation of planetesimals probably requires periods (and/or regions) of 
quiescence such as dead zones or pressure bumps (see Kretke & Lin 2007; Brauer et al. 
2008b). These locally confined quiescent environments could decrease both the relative 
and radial velocities of particles, thus allowing particles possibly to either grow to larger 
sizes or to accumulate until gravitational instabilities become important (see Johansen 
et al. 2007; Lyra et al. 2009). These effects are not yet included in this model and will 
be the subject of future research. 

Assuming that the proposed mechanism of dust retention is effective in protoplan- 
etary disks, we can constrain the critical fragmentation velocity to be smaller than 10 
m/s since velocities > 10 m/s cannot explain the observed lifetimes of disks even in the 
case of small a. 
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Figure 5.2: Dust-to-gas ratio (solid line, linear scale) and 0.55 /im optical depth 
at 10 AU (dashed line, logarithmic scale) as function of time for model 1 (top), 
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CHAPTER 

Testing the theory of grain 

growth and fragmentation by 

millimeter observations of 

protoplanetary disks 

Published as Birnstiel, Ricci, Trotta, DuUemond, Natta, Testi, Dominik, 
Henning, Ormel, & Zsom 201 Oh, A&A, 503, L5 

Abstract 



Observations at sub-millimeter and mm wavelengths will in the near 
future be able to resolve the radial dependence of the mm spectral slope 
in circumstellar disks with a resolution of around a few AU at the distance 
of the closest star-forming regions. We aim to constrain physical models of 
grain growth and fragmentation by a large sample of (sub-)mm observations 
of disks around pre-main sequence stars in the Taurus- Auriga and Ophiuchus 
star-forming regions. State-of-the-art coagulation/fragmentation and disk- 
structure codes are coupled to produce steady-state grain size distributions 
and to predict the spectral slopes at (sub-)mm wavelengths. This work 
presents the first calculations predicting the mm spectral slope based on a 
physical model of grain growth. Our models can quite naturally reproduce 
the observed mm-slopes, but a simultaneous match to the observed range of 
flux levels can only be reached by a reduction of the dust mass by a factor 
of a few up to about 30 while keeping the gas mass of the disk the same. 
This dust reduction can either be caused by radial drift at a reduced rate 
or during an earlier evolutionary time (otherwise the predicted fluxes would 
become too low) or due to efficient conversion of dust into larger, unseen 
bodies. 
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6. Testing the theory of grain growth and fragmentation 
6.1 Introduction 

Circumstellar disks play a fundamental role in the formation of stars as most of the 
stellar material is believed to be transported through the disk before being accreted 
onto the star (Lynden-Bell & Pringle 1974). At the same time circumstellar disks are 
thought to be the birth places of planets. Understanding the physics of circumstellar 
disks is therefore the key to some of the most active fields of astrophysical research today. 

However, observing these disks in order to learn about the physical processes taking 
place in their interior is a challenging task. Strom et al. (1989) and Beckwith et al. 
(1990) were the first to use observations at mm-wavelengths to confirm that many of the 
observed pre-main sequence (PMS) stars showed excess radiation above the spectrum of 
a T Tauri star. While these single dish observations provided valuable insight into dust 
masses (because mm observations probe not only the thin surface layers, but the bulk of 
the dust mass in the mid-plane), recent developments in the field of (sub-)mm interfer- 
ometry allow one to constrain models of disk structure and evolution of protoplanetary 
disks by fitting parametric models to the observed radial profiles (e.g., Andrews et al. 
2009; Isella et al. 2009). Spatially resolving the disks is important because it ensures 
that low millimeter spectral slopes are not just an artifact of high optical depth. 

Today, mm spectral slopes are known for quite a number of disks, and spatially 
resolved observations indicate that the low values measured in these samples are related 
to grain growth (e.g., Testi et al. 2003; Natta et al. 2004; Rodmann et al. 2006). Grains 
are believed to collide and stick together by van der Waals forces, thus forming larger 
and larger aggregates (Dominik &: Tielens 1997; Poppe et al. 2000; Blum & Wurm 2008). 
Due to this loose binding, collisions with velocities in excess of a few m s^^ may lead to 
fragmentation of the aggregates. 

Larger samples of radially resolved mm spectral slopes are expected in the near 
future, but so far no study interpreted mm observations using simulated grain size dis- 
tributions. Instead simple parametric power-law distributions were used with an upper 
size cut-off. In this work, we use a state of the art dust grain evolution code (similar 
to Brauer et al. 2008a; Birnstiel et al. 2010a) to derive steady-state grain distributions 
where grain growth and fragmentation effects balance each other. We self-consistently 
solve for the grain size distributions and the disk structure to predict fluxes at mm wave- 
lengths and the radial dependence of the mm spectral index. Comparing these results 
with observed values in the Taurus and Ophiuchus star-forming regions allows us to test 
the predictions of the theory of grain growth/fragmentation and to infer constraints on 
grain properties such as the critical collision velocity and the distribution of fragments 
produced in collision events. 

Grains orbiting at the Keplerian velocity in a laminar gas disk feel a constant head 
wind (caused by the gas rotating slightly sub-keplerian) , which forces them to spiral 
inwards (Weidenschilling 1977). If this drag is as efficient as the laminar theory predicts, 
all dust particles that are necessary to explain the observed spectral indices would quickly 
be removed (see Brauer et al. 2007). We therefore assume that the radial drift is halted 
by an unknown mechanism. Under this assumption, we find that the low values of the 
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6.2 Model description 



Table 6.1: Parameters of the model grid: Mjisk is the total disk mass, a^ is the 
turbulence parameter, Ui is the critical collision velocity, /vac is the grain volume 
fraction of vacuum, and ^ is the index of the distribution of fragments (see Eq. 6.3). 
The parameters of the fiducial model are highlighted in bold face. 
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mm spectral index can be explained by the theory. We show that in order to explain 
the observed flux levels, the amount of observable^ dust needs to be reduced by either 
reducing the dust-to-gas ratio (perhaps by radial drift at a intermediate efficiency or 
during an earlier evolutionary epoch) or by dust particle growth beyond centimeter 
sizes. 



6.2 Model description 
6.2.1 Disk model 

We consider disks around a PMS star with a mass of 0.5 M©, bolometric luminosity of 
0.9 Lq and effective temperature of 4000 K, at a distance of 140 pc, which are typical 
values for the sample of low-mass PMS stars studied in Ricci et al. (2010, hereafter RIO). 
To derive the disk structure we adopted a modified version of the two-layer models of 
passively irradiated flared disks developed by Dullemond et al. (2001) (following the 
schematization by Chiang & Goldreich 1997), in which we have relaxed the common 
assumption that dust grain properties are constant throughout the disk. For the disk 
surface density we adopted the self-similar solution for a viscous disk (see Lynden-Bell 
& Pringle 1974) with parameters lying in the ranges observationally constrained by 
Andrews et al. (2009). The surface density gradient 7 and the characteristic radius R^ 
(for the definitions, see Hartmann et al. 1998) are assumed to be 7 = 1 and i?c = 60 AU, 
respectively. Throughout this work we assume a constant dust-to-gas mass ratio of 1%. 



^by "visible" or "observable" dust we mean the dust particles which are responsible for most 
of the thermal continuum emission at (sub-)mm wavelengths, which are typically smaller than a 
few centimeter in radius. 
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6. Testing the theory of grain growth and fragmentation 
6.2.2 Dust model 

We use a coagulation/fragmentation code as described in Brauer et al. (2008a) and 
Birnstiel et al. (2010a) to simulate the growth of dust particles. Particles grow through 
mutual collisions (induced by Brownian motion and by turbulence, see Ormel & Cuzzi 
2007) and subsequent sticking by van der Waals forces. We assume the dust particles to 
be spheres of internal density ps and vary p^ to account for porosity effects. However, 
we do not employ a dynamic porosity model (see Ormel et al. 2007; Zsom & Dullemond 
2008). 

With increasing collision velocity Am, the probability of sticking decreases and frag- 
mentation events start to become important. Here, we use the fragmentation probability 



Pf 



if Au < Uf — 5u 

1 ifAu>Uf , (6.1) 

ou 

where Uf is the collision velocity above which particles are assumed to fragment, and 
5u is the transition width between coagulation and fragmentation (taken to be 0.2 Uf). 
Recent studies of collision experiments (Giittler et al. 2010) and numerical simulations 
(Zsom et al. 2010) indicate that there is also a regime in which particles may bounce. 
However, because this topic is still not well understood, we will omit these effects here. 

Radial drift is an as yet unsolved problem (Brauer et al. 2007; Birnstiel et al. 2009). 
Still there are several effects such as spiral wave structure (e.g., Cossins et al. 2009), 
density sinks (e.g., Brauer et al. 2008b), or zonal flows (e.g., Johansen et al. 2009) which 
may reduce the effectiveness of radial drift. We assume that radial drift is ineffective 
because we focus on the question whether observations can be explained through the 
physics of grain growth and fragmentation. The question to answer in this case is not 
how to retain particles at these radii, but rather how to create them there in the first 
place. 

To investigate this problem, we simulate the physics of particle growth and frag- 
mentation until a steady state between both processes develops. Because the relative 
velocities for particles typically increase with grain radius, we can relate the fragmen- 
tation velocity to a certain grain size (which defines the "fragmentation barrier", see 
Birnstiel et al. 2009) 



2S, u^ 



' (6.2) 



vratPs ci 

above which particles fragment (with Sg, at and Cg being the gas surface density, the 
turbulence parameter, and the sound speed, respectively). Uf and at are assumed to be 
radially constant with at values within a range expected from theoretical (see Johansen & 
Klahr 2005; Dzyurkevich et al. 2010) and observational works (see Andrews et al. 2009). 
Grains which reach Omax will experience high velocity collisions, causing them to be 
eroded or even completely fragmented. The resulting fragments can again contribute to 
growth processes at smaller sizes, and the grain size distribution will at some point reach 
a steady state where gain and loss terms caused by coagulation and by fragmentation 
cancel out at all sizes. 
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Particles will need a certain time to grow to the fragmentation barrier. The time to 
reach the steady state will therefore be several of these growth time scales. Depending 
on the distance to the central star, the steady state is typically reached after a few 
thousand years at 1 AU up to about 1 Myr at 100 AU. The mean ages of the sources in 
our sample are ~ 2 Myr and k 0.5 — 1 Myr for the Taurus and Ophiuchus PMS stars, 
respectively. Since radii around 40 — 80 AU dominate the observed emission at (sub-)mm 
wavelengths, we expect most of the samples to be in or at least close to a steady state. 

If the highest collision velocity that turbulent motion induces (depending on at and 
Cs) is lower than the critical collision velocity Uf, then (at least some) particles do not 
fragment (i.e. the break through the fragmentation barrier) and a steady state is never 
reached. Owing to this scenario, some of the possible combinations of the parameter 
values (see Table 6.1) do not reach a steady state and are therefore not included in the 
results. 

The shape of the steady-state grain size distributions is influenced mainly by five 
parameters: the previously mentioned at, Uf, Sg, the temperature T (through the sound 
speed Cs), and by the prescription of fragmentation. In our models, we assume the 
distribution of fragments to follow a power-law number density distribution, 

n{m)ccm^^ , (6-3) 

with an upper end at nif. We consider fragmentation and cratering, as described in 
Birnstiel et al. (2010a). Recent experiments suggest ^ values between 1.07 and 1.37 (see 
Giittler et al. 2010). In this work, we consider ^ values between 1.0 and 1.8. 

To calculate the dust opacity of a given grain size distribution we adopted the same 
dust grain model as in RIO, i.e. porous composite spherical grains made of astronomical 
silicates, carbonaceous materials and water ices (see RIO for the references to the optical 
constants). The ratio between the fractional abundance of each species comes from 
Semenov et al. 2003, and models with three different porosities have been considered 
in this chapter (see Table 6.1). We used the Bruggeman mixing theory to combine 
the refractive indices of the different materials and to calculate the dust opacity of 
the composite grains. The opacity induces probably the largest uncertainties in our 
calculations because grain composition, grain structure and temperature effects may 
lead to largely different opacities (see, for example Henning &: Stognienko 1996). 

6.2.3 Comparison to observations 

We compare the (sub-)mm SED generated by our models with observational data of 
RIO and Ricci et al. (in prep), more specifically the flux at 1 mm (Fimm) and the 
spectral index between 1 and 3 mm (F(A)ocA^"i"^™™). The samples considered include 
all class II disks in the Taurus-Auriga and p-Oph star-forming regions for which both 
the central PMS star and the disk are observationally well characterized through optical- 
NIR spectroscopy/photometry and (sub-)mm photometry/interferometry. To calculate 
the dust opacity as a function of wavelengths and radius, and the temperature in the 
disk mid-plane, we iterated the two-layer disk model (keeping the profile of Eg constant 
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6. Testing the theory of grain growth and fragmentation 



in time) with the dust model described above until convergence is reached. Once the 
physical structure of the disk is determined, the two-layer disk models return the disk 
SED, which can be compared with the observations. 

The influence of the different parameters on the calculated ai-3mm values can mostly 
be understood by a simple model for a dust distribution, as used in RIO (cf. Fig. 3 
in RIO): for maximum grain sizes much smaller than the observed wavelengths, the 
spectral index of the dust opacity /3i_3mm (k(A)ocA^'^I"^™™) is constant, while it decreases 
for maximum particle sizes larger than sub-mm. In between (at a few tenth of a mm), 
there is a peak which is caused by an increased opacity of grains with sizes similar to the 
observed wavelength. The relation between ai_3mm and /3i-3mm depends on the emitting 
spectrum and the optical depth. For a completely optically thin disk in the Rayleigh- 
Jeans regime /3i_3mm = ai-3mm — 2. However, if the emitted spectrum deviates from the 
Rayleigh-Jeans limit, then Pi.^ram ^ cti-Smm — 2. In our models, ai_3nim — /3i-3mm turns 
out to be typically between 1.4 and 1.7 if Omax is outside of the peak of opacity. 



6.3 Results 

6.3.1 Sub-mra fluxes and spectral indices 

For all possible combinations of the parameters shown in Table 6.1, we solved for the 
steady-state grain size distributions and derived the ai-smm and -Fimm values. As noted 
before, some of the models do not result in a steady state and are therefore not shown 
here. 

The top left panel of Fig. 6.1 shows the influence of the turbulence parameter at- 
According to Eq. 6.2, the maximum grain size increases if at decreases. Depending on 
where amax lies with respect to the opacity peak (see RIO, Fig. 3), ai-3mm can increase 
or decrease with increasing at. In the simulations presented here, Omax is typically so 
large that increasing at predicts steeper spectral slopes. 

flmax is more sensitive to Uf (cf. top right panel in Fig. 6.1): the maximum grain size 
flmax is proportional to uf, therefore a change of U{ by a factor of about 3 significantly 
changes ai_3mm by increasing the grain size by about one order of magnitude. However 
many models with a fragmentation velocity of 10 m/s never reach a steady state. It is 
therefore not possible to explain lower ai_3mm values by a further increase of Uf alone. 

The influence of Mjisk on Fimm and ai-3mm is twofold. Firstly, a decrease in Mdisk 
(assuming a constant dust-to-gas ratio and a fixed shape of the disk surface density, i.e. 
not varying Re and 7, see Section 6.2.1) reduces the amount of emitting dust and thus 
-^imm- Secondly, such a reduction in gas mass also reduces Omax (Eq. 6.2), which tends 
to increase ai_3nim- This combined trend is seen in Fig. 6.1. Hence, in order to explain 
faint sources with low ai-smmj the amount of emitting dust has to be reduced while the 
disk gas mass stays large. This effect could be achieved in two ways: the amount of 
dust could be reduced by radial drift at a reduced rate (full radial drift would quickly 
remove all mm-sized grains, see Brauer et al. 2007) or only the "visible" amount of dust 
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Figure 6.1: Influence of the parameters at (top left), fragmentation velocity (top 
right), disk mass (bottom left) and grain porosity (bottom right) on the observed 
fluxes and spectral indices. The black circle denotes the flducial model whose 
parameters are given in Table 6.1. The grey area represents the region in which 
the observed sources lie (see Fig. 6.2). 



is reduced if some of the dust is already contained in larger bodies. This latter case 
is predicted by our non-steady-state distribution models and will be discussed in more 
detail in a forthcoming paper. 

In general, lower values of ^ translate to shallower grain-size distribution, which 
results in lower values of /3i_3mm (see Draine 2006). The lower right panel in Fig. 6.1 
does not seem to indicate a strong dependence on ^, but lower values of ^ (around 1) 
seem to be closer to the observations especially at high fluxes. 

Figure 6.2 shows the areas which are covered by our sets of simulations for different 
porosities in comparison to the observational samples. It can be seen that only the 
brightest sources are covered by the simulations. The trend of larger ai_3mm for a larger 
vacuum fraction seems to be in contradiction with Eq. 6.2, because smaller grain volume 
density leads to larger Omax- However in this case, the opacity is much more affected by 
changing the grain structure: reducing the grains' vacuum fraction increases the spectral 
index at sub-mm, while it is reduced for longer wavelengths. Therefore opacity effects 
outweigh the smaller changes in Omax- A more thorough analysis of opacity effects is 
beyond the scope of this chapter, but it seems implausible that the large spread in the 
observations could be explained by different kinds of grains alone (see Draine 2006). 
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Figure 6.2: Observed fluxes at mm-wavelengths of the Taurus (black dots) and the 
Ophiuchus (blue dots) star-forming regions (see RIO and Ricci et al. 2010, in prep) 
and the areas covered by the simulation results for different vacuum fractions of 
the grains (varying all other parameters according to Table 6.1). 



6.3.2 Radial profiles of the dust opacity index 

The presented models also compute ai_3mm as function of radius. Prom the point of view 
of a comparison with the observations, this is somewhat premature because observational 
methods are not yet able to provide reliable radial profiles of ai-3mm (e.g., Isella et al. 
2010, Banzatti et al., in prep.). Yet, the predicted radial dependence of /3i-3mm (shown 
in Fig. 6.3) agrees with the observations so far. It can be seen that the shape of most 
models looks similar, slightly increasing from /3i_3inin-values around 0.5 at 10 AU up to 
around 1.5 at 100 AU. The reason for this is that Omax depends on the ratio of surface 
density over temperature. Under typical assumptions, Omax will decrease with radius. 
An upper grain size, which is decreasing with radius and stays outside the peak in the 
opacity, results in f3i.3rara increasing with radius (cf. Fig. 3 in RIO). If the radially 
decreasing upper grain size Omax reaches sizes just below mm, then the peak in opacity 
will produce also a peak in the radial profile of /3i-3mm (the size of which depends much 
on the assumed opacity), which can be seen in Fig. 6.3. Thus, even though Omax is 
monotone in radius, /3i_3mm does not need to be monotone. 
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Figure 6.3: Predicted profiles of the dust opacity index at mm-wavelengths for 
different variations of the fiducial model. The colors correspond to the parameters 
shown in Fig. 6.1. 



6.4 Discussion and conclusions 

We present the first in-depth comparison of simulated grain size distributions and ob- 
served mm spectral indices of YSOs in the Taurus and the Ophiuchus star-forming 
regions. Additionally we present the first predictions of the radial profile of the dust 
opacity index at mm wavelength that are consistent with the limits set by Isella et al. 
(2010). 

Low values of the observed mm-slopes are quite naturally reproduced by our mod- 
els, favoring low values of S, and at as well as fragmentation threshold velocities above 
1 m s^^. However, a simultaneous match to the observed range of fiux levels requires a 
reduction of the dust mass by a factor of a few up to about 30. This over- prediction of 
fiuxes cannot be fixed by simply reducing the disk mass because the predicted ai_3mm 
would be too large for smaller disk masses. Opacities induce a large uncertainty in the 
flux levels. However, considering the results of Draine (2006), it seems implausible that 
the large spread in observed fluxes for different disks with similar ai_3mm (which is prob- 
ably even larger as very faint disks are not contained in the sample) can be explained 
by different grain mineralogy alone. 

The aforementioned reduction of observable dust could be due to radial drift at a 
reduced rate or during an earlier epoch (drift has been artificially suppressed in this 
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work in order to explain the low values of ai-smm by > 1 mm sized grains). Another 
possible explanation is grain growth to even larger sizes, as these bodies have a small 
opacity coefficient per unit mass. 

Finally, a different dependence between ai_3mm and the observed flux Fimm might 
also originate from disk surface densities profiles that differ from what we have assumed 
in this work. This possibility, as well as the effect of a different dust composition, will 
be considered in a future work. 
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CHAPTER 




Summary and outlook 



In this thesis, we investigated how dust in circumstellar disks evolves due to radial 
transport within the gaseous disk and due to grain growth and fragmentation. The 
physics of the underlying gas disk is modeled as well using a viscous accretion disk 
model including irradiation by the central star and viscous heating. 

The first two chapters of this thesis focused on the various physical aspects of the 
evolution of dust and the gas in these disks using a newly developed code which can 
integrate the gas disk evolution and at the same time the coupled equations of advection, 
diffusion and the Smoluchowski equation including fragmentation of grains. Additionally, 
the solver of the Smoluchowski equation can in the future be extended to also solve for 
the mean porosity of each grain size, as derived and demonstrated by Okuzumi et al. 
(2009). 

Our results show that the radial inward drift and the fact that larger grains are less 
effectively mixed and carried outwards causes a strong radial dependence of the dust- 
to-gas ratio. Only if grain fragmentation effectively grinds down larger grains, the dust 
can be relatively equally be distributed with the gas. Observations in the near future 
(e.g., ALMA) should be sensitive enough to detect the different sizes of the gas and the 
dust disk (Panic & Birnstiel, in prep.). 

We also investigated the two main barriers preventing particles to accumulate to 
planetesimals: radial drift and fragmentation. Our results show that there are several 
mechanisms which are able slow down radial drift to an extend where pure sticking would 
allow particles to overcome the radial drift barrier. However, fragmentation seems to be 
a much stronger obstacle, especially in the inner hot parts of the disk. A very low degree 
of turbulence is needed to keep the relative velocities of the particles low enough to avoid 
shattering. Even in this case, radial and azimuthal relative velocities are typically still 
large enough to fragment the aggregates. Therefore, we propose that both a low degree 
of turbulence and a low radial pressure gradient are necessary to allow further growth 
by sticking collisions. 
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Dead zones seem to be such an environment. Dzyurkevich et al. (2010) found that 
indeed a pressure bump can exist at the inner boundary of a dead zone. Brauer et al. 
(2008b) showed that particles can break through the growth barrier at a pressure bump. 
Preliminary simulations of such a pressure bump, caused by a radial variation in the 
amount of turbulence in an evolving gas disk, show that radially drifting grains accu- 
mulate in the center of the pressure peak where the physical conditions (low turbulence, 
vanishing radial and azimuthal relative velocities) are such that larger grains grow fast 
enough to overcome the growth barrier. A snapshot of such a simulation is shown in 
Fig. 7.1. Future work in this direction has to take into account more detailed profiles 
of the relative velocities and of the gas surface density, such as produced by MHD sim- 
ulations Dzyurkevich et al. (2010). The code presented in this thesis can possibly be 
extended to a layered disk scenario. These simulations could give some first results about 
the enrichment in dust and the dust size evolution within a dead zone. 



time — 1.19 x 10^ years 




-10 



Figure 7.1: An underlying change in the turbulence parameter at (from 10^"^ inside 
30 AU to 10^^) induces a surface density and thus pressure maximum. Dust 
accumulates in this pressure bump and - because radial and azimuthal velocities 
vanish there - dust can continue to grow to larger sizes. The white line denotes 
the grain size which corresponds to a Stokes number of unity and is proportional 
to the gas surface density. 



A second step in this direction would be to use the code to solve not the radial but 
instead the vertical transport of dust caused by settling and turbulent mixing. This could 
be coupled to the already available vertical temperature and density profile of the disk 
to investigate the observable dust distribution in the upper layers of disks. Together 
with a dust-dependent recipe for the MRI activity (available from Neal Turner), this 
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would allow us to study the interplay between MRI activity (which in turn depends on 
the amount of small dust) and the dust evolution (which depends on the MRI activity 
as source of turbulent velocities). 

Another result of this thesis was that if radial drift is ineffective or if particles stay 
small enough (due to the fragmentation), then grain growth and fragmentation balance 
each other. This balance determines the size distribution of dust in circumstellar disks. 
Therefore, in Chapter 4, we tried to analytically understand the dust distribution in 
such a case. For this purpose, we generalized previous analytical results for collisional 
cascades to include both processes, coagulation and fragmentation at the same time. We 
could derive power-law distributions for three different regimes, depending on which of 
these effects dominates the shape of the distribution. 

For simple collision kernels, we found very good agreement between the analytical 
results and the numerically simulated steady-state distributions. For the more compli- 
cated collision kernel which represents the condition in a turbulent circumstellar disk, 
our analytical results helped us to understand the general shape of the resulting distri- 
butions. We were able to find a recipe which can reproduce the distributions we derived 
numerically. This recipe can readily be used for further modeling, such as grain surface 
chemistry (Vasyunin et al., in prep.), disk heating-cooling and FUV photoevaporation, 
radiative transfer modeling of disks (Mulders &: Birnstiel, in prep.; Sauter &: Birnstiel, 
in prep.) and many more. 

As mentioned before, the dust distribution is an important ingredient to many models 
of circumstellar disks and is thus a key to understanding the astrophysics of circumstellar 
disk and of planet formation. Chapters 5 and 6 are an example of this: in Chapter 5, 
we showed how grain fragmentation and cratering could solve the riddle why disks are 
observed to be dusty for several million years, while theory predicts the dust to grow 
and to disappear quickly due to radial drift: if fragmentation is happening at an collision 
velocity of around 1 m s^^, as suggested by experiments (Blum &: Wurm 2008), then 
the maximum grain size which can exist is small enough to be unaffected by radial drift. 
Particles larger than this size experience larger impact velocities and are quickly ground 
to smaller sizes. 

Another link to observations was established in Chapter 6: we argued that some 
mechanism must exist to stop radial drift. Otherwise, all larger grains in the outer disk 
would disappear quickly, as shown by Brauer et al. (2007). This assumption allows us to 
investigate the question whether grains can grow at these radii to the observed sizes on 
a reasonable timescale. We found that the grain size distribution reaches a steady-state 
after about a million years, even at 100 AU. The largest grain size is therefore not limited 
by the growth time scale, but by fragmentation. We calculated disk models for a large 
set of models to determine the flux and the spectral index predicted by the theory and 
compared it to a sample of observations by Ricci et al. (2010). We found that the models 
indeed reproduce the low spectral indices, however the fluxes are only reproduced for the 
brightest objects. This indicates a reduced dust-to-gas ratio (or decreased dust opacity) 
since increasing both the dust and gas mass would increase the spectral index. Future 
work in this direction will focus on the impact of using other opacity models and on 
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time evolving simulations of grain growth and fragmentation. Preliminary results in this 
direction show that disk at earlier evolution times produce lower flux values, as shown 
in Fig. 7.2. 
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Figure 7.2: Preliminary results: comparison of fluxes and spectral indices of ob- 
servational data (black and blue dots for Taurus and Ophiuchus, respectively), 
steady-state solutions from Chapter 6 (the area inside the red polygon) and time 
evolving simulations (green crosses). The snapshots are taken at 10^ years of evo- 
lution. Higher-mass disks are already inside of the steady-state area, while lower 
mass disks reproduce some of the lower-fiux values. 

It has now been mentioned several times, that planet (-esimal) formation heavily 
depends on this initial phase of grain growth. We will therefore use the code presented 
in this thesis also to study the inward mass flux of dust and the maximum size of the 
drifting particles (Birnstiel & Klahr, in prep.). The reasons for this are the following: 
if planetesimals are to be formed by gravoturbulent mechanisms (see Johansen et al. 
2007), then the initial conditions for these scenarios need to be met. The conditions 
are firstly that grains need to have the right sizes (for example about decimetre sized 
particles for the scenario of Johansen et al. 2007) and secondly that there needs to be 
enough of these particles (i.e., a strongly enhanced dust-to-gas ratio). Our model can 
calculate self-consistently, which particles are produced under which condition and at the 
same time the radial drift and possible accumulation of the dust. This drift rate and the 
particle sizes are also important for other alternative planet formation scenarios, such 
as dust capturing in vortices (Barge &: Sommeria 1995; Klahr & Bodenheimer 2006). 
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